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A NUMERICAL STUDY OF ELECTROMAGNETIC SCATTERING 
FROM OCEAN-LIKE SURFACES 

R„ R 0 Lentz 

The Ohio State University ElectroScience Laboratory 

CHAPTER I 
INTRODUCTION 


The scattering of electromagnetic waves from the ocean surface 
has been of great interest for some time. In this work the scat- 
tering from one dimensional sea-like random surfaces is ex- 
amined by a variety of computational methods, with a view to 
establishing what practical limitations must be satisfied on such 
surface parameters as radius of curvature, mean squared height, etc., 
in order that the statistical properties of the scattered radiation 
may be calculated with reasonable accuracy. The results of the com- 
putations are then used to discuss the applicability of the several 
theoretical models for sea-surface scattering (geometrical optics, 
physical optics, perturbation theory and the composite model) and the 
prospect for direct calculation of the scattered fields from the 
actual sea surface. 

During the past few years, theoretical and experimental work here 
and abroad (Refs. [l]-[7]) has led to an understanding of the mech- 
anisms responsible for scattering and emission of microwaves by the 
ocean. For off-normal backscatter, the "Bragg-scatter" from capil- 
lary and short wavelength components of the ocean surface, which 
can be calculated by perturbation theory, has explained the angular 


and polarization dependence of the microwave radar return. When 
combined with the known height spectrum (Ref. [8]) of the ocean 
surface, it explains the weak dependence of backscatter on electro- 
magnetic wavelength and wind velocity. Near the specular direction, 
i.e., near normal incidence for backscatter, the scattering is con- 
trolled by the slope distribution of the large scale structure of 
the surface. This part of the scattering is calculated by geometrical 
optics, and explains the dependence of the emissivity of the surface 
on wind velocity. 

Nevertheless, the many assumptions required in finding the 
scattered fields by the perturbation or geometrical optics approxi- 
mations, particularly assumptions about the Gaussian character of the 
surface height statistics, and the applicability of the theoretical 
approximations to the actual sea surface, have led to considerable 
discussion about the validity of the various theoretical solutions 
(Ref. [9]). Since straightforward verification by measurement is 
not practical, partly because of difficulty in the measurement process 
itself and partly because of the difficulty in specifying exactly what 
the surface was when the measurement was being made, it is desirable to 
have a direct method for calculating the scattering from a specific 
realization of the ocean surface. Direct calculations will allow a 
realistic assessment of the validity of the various theories, without 
any assumptions about the statistical properties of the surface. If 
a statistical average of the scattered fields over an ensemble of 
surface representations is required, it can be obtained (albeit at 



some cost) by a direct summation of the scattered fields from the 
individual surface representations. 

The specific surfaces considered here are cylindrical perfectly 
conducting surfaces as shorn in Fig. 1. The surface generators are 



Fig. 1.— The scattering surface. 


parallel to the z axis, and the surface elevation is specified by 
y = H{X). The incident field is a plane wave whose direction of 
propagation lies in the x 3 y plane and makes an angle of THI with 
the positive x axis., while the observation direction makes an angle 
of THS with the positive x axis. Time dependence is assumed to be 
e JtI)i: and has been suppressed throughout. All distances are measured 
in centimeters. 

Three different methods for calculating the fields from such 
a surface are developed here. Although the details are discussed 
later it is desirable to outline each technique at this time. 


The first approximate method is the geometrical optics tech- 
nique (G.O.). For a given surface, and given scattering and in- 
cidence angles, the program locates the specular points on the 
surface (points where the local incidence angle equals the local 
scattering angle) and evaluates the radius of curvature at each 
specular point. The scattered far field is then found by summing 
the contribution from each of the specular points, including an 
extra 90° phase shift for the fields scattered from concave up 
portions of the surface. Shadowing of one section of the surface by 
another section may be taken into account. 

The next approximation is the physical optics (P.0.) technique. 
For a given surface the scattered field is computed by integrating 
over the approximate surface current 

O ) J s = 2n x H 1 

where n is the outward normal to the surface and H 1 is the incident 
magnetic field. Shadowing is always taken into account, as this is 
implicit in the physical optics formulation. 

The last method developed here is based on a point matching 
solution to the integral equation satisfied by the true surface 
current J* s . The scattered fields are then found by integrating 
over the surface currents. Test cases (e.g., the wedge 
problem) have shown this method to be by far the most accurate; 
hence it is used as a standard to which all others are compared. 
However, because of computer storage limitations, this program can 
not handle surfaces whose arclengths are greater than ^60 electrical 



v/avelengths , v/hereas the G.O. and P.0, programs can, in principle, 
handle surfaces of any length provided sufficient computer time is 
available. 

In order to avoid edge effects, tapering of the incident field 
is necessary in the integral equation solutions. The same tapering 
has been applied in both the G.O. and P.0, solutions so that they 
can be directly compared to the exact fields. The tapering applied 
here is illustrated in Fig. 14 of Chapter IV. 

In the succeeding chapters each of these methods will be 
described in detail. By comparing the results for a series of 
test surfaces, the limitations of each method are established. 


CHAPTER II 


THE GEOMETRICAL OPTICS METHOD 

The first approach to examining the scattering from a one 
dimensional rough surface is the geometrical optics method. By 
this is meant that the scattered field is computed by finding 
the specular points on the surface, and associating with each such 
point a scattered field amplitude and phase which depend on the 
geometrical properties of the surface at the specular point. 

A. Geometrical Optics 

Conservation of energy flux along a ray path will provide us 
with the geometrical optics field strengths (Ref. [10]). 

Consider the two dimensional ray tube shown in Fig. 2. If u Q 
is the field strength at some reference point at a distance p from 
the caustic and u is the field strength at distance p + z from the 
caustic, then the conservation of energy in the ray tube requires 



Fig. 2. — Ray tube geometry. 


(2) u 2 p de = u 2 (p + a) do 

so that one may write 


(3) 


uW " 


The factor e”^ <5 ' 9 with A 


the electrical wavelength and 


(4) k = 27r/A e 

accounts for the phase shift between p and p+A. Equation (3) fails 
at £ equal to -p. This location (at the confluence of the rays) is 
termed a caustic. Kay and Keller (Ref. [11]) have demonstrated that 
at points beyond the caustic (& less than -p) Eq. (3) is still valid 
if a phase shift of +90° is introduced. 

To use geometrical optics it is necessary to find all points 
on the scattering body at which the law of reflection is satisfied 
locally for the particular set of THI and TH5 under consideration. 

Once these points are located Eq. (3) is used to calculate the scat- 
tered field. Figure 3 shows the geometry for the calculation of 
the scattered field from one such specular point. By the law of 
reflection 3 the local incidence and scattering angles are equal and 
are marked ANG in the figure. The distances marked r c and p are the 
radius of curvature and the distance from the specular point to the 
optical image of the source (i.e. 9 the caustic distance) respectively. 
The distance p is given by a cylindrical mirror formula as 
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Fig. 3.— Specular point geometry. 


(5) 


1 _ 2 ,1 
p' |r c | cos (ANG) * 


In the cases considered here the distance to the line source, z q9 
will be assumed to be infinite, hence 

| r _ | cos (ANG) 

(6) p = 5 . 


If the specular point is taken as the reference position then Eq. (3) 
gives u $ , the scattered field at the observation position 


(7) 


u = R u. . 
s 1 \l 


-jk£ 


p+Z 


- R u^ ,fp e~^ z /yJT for z » p (far field) 
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where is the incident field evaluated at the specular point and 
R is a reflection coefficient. If the electric field is parallel 
to the surface generators (T.M. case) and u^ is taken as the inci- 
dent electric field,, then u s is taken as the scattered electric 
field with R = -1 . If the magnetic field is parallel to the surface 
generators (T.E. case) and u^ is taken as the incident magnetic fields 
then u s is the scattered magnetic field and R = +1 . For dielectric 
scatterers the corresponding Fresnel reflection coefficients are to 
be used for R. This makes the geometrical optics program the easiest 
to convert from perfectly conducting bodies to penetrable bodies. 

Up to this point the scattering surface has been assumed to 
be concave down at the specular point. If the body is concave up 
at the specular point then the caustic position is above the surface 
instead of below, the scattered rays pass through the caustic on 
the wa y to the observation point if the observer is in the far field, 
and thus a phase shift of +90 degrees must be introduced. The 
distant scattered fields may then finally be written 


( 8 ) 



Specular 

Point 


r | cos (ANG) 

V 

2 


e" jk * 

A 


for the T.M. case and 

h|U) = 4 


r | cos (ANG) 


e “jk.e 


( 9 ) 


Specular 

Point 


e 


for the T.E. case, where e is +1 if the surface is concave down at 
the specular point and +j if the surface is concave up at the specular 
point. 

On an actual surface there may be several specular points con- 
tributing to the total scattered field, so it is important to pre- 
serve the phase relationships among them. Phase reference is taken 
at the origin, and an incident wave of unit amplitude is assumed, 
i .e. , 


(10) E’ = e' J ’ k ‘ R (T.M. case) 

(11 ) H’ = e' jk ' R (T.E. case) 

where 

(12) k • R = |2_(-x cos (THI) - H(X) sin (THI)). 

A e 

With the aid of the geometry shown in Fig. 4, the scattered far 
field is found from Eqs. (8) and (9), with i = where 

(13) = -R * D s = - x cos (THS) - H(x) sin (THS), 
and 

(14) D s - cos(THS)x + sin (THS) y 


10 




The total scattered field in the THS direction is the sum of the 
fields scattered by each of the specular points. The numerical 
values of the scattered fields as calculated by the programs of 
Appendix A, and plotted in the various figures of Chapter V are 
denoted by e| and H^ 9 and have been normalized with respect to 
the actual fields e|(^),h|(^) by 

E z<V 

>- 

H^) . 

> 

It is clear that Eqs. (15) and (17) fail if the radius of 
curvature is infinite at the specular point. This is because the 
source was assumed at infinity, i.e., °°. If $, were to be 

held finite then from Eq. (5) 

(19) lim P = £ 

o 

r -+■ oo 

c 

and the singularity in Eqs. (15) and (17) would not occur. In ad- 
dition to the singularities caused by an infinitely distant source, 
there are a number of other shortcomings of the G.O. approximation. 
Among them are: a failure to account for wedge diffraction effects 

(radius of curvature goes to zero), a failure to account for dif- 
fraction from shadow boundaries into shadowed regions (Ref. [12]), 
a failure to properly predict the scattered fields if the surface 



r s ^ 
E, 



z j 

jk£, 

> = /sy e 1 - 

(18) -< 

H z 

J 


features subtend only a few Fresnel zones (Ref. [13]),and finally 
a failure to predict any scattered field if no specular point exists 
on the body. 

Implicit in the geometrical optics technique is the concept of 
shadowing, that is, a specular point cannot contribute to the scat- 
tered field unless it can be seen by both the source and the observer. 
The program developed here can account for shadowing of this type. 

B. Discussion of the Geometrical Optics Program 

For geometrical optics calculations the first order of business 
is the location of the specular points. Figure 5 shows the geometry. 


y 



The surface height profile is described by H(X) and the regions under 
investigation lies between XSTRT (X START) and XSTOP. THI and THS 
have already been defined; THN (THETA of the NORMAL) is the angle 

A 

between the normal (n) to the surface and the positive x axis. 

Clearly 
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( 20 ) 


THN(x) = ir/2 + Tan " 1 (dH(x)/dx). 

The law of reflection gives (x,H(X)) as a specular point when 

(21) THS - THN(X) = THN(X) - THI 
i .e. 3 

(22) (THS + THI )/2 = THN(X) . 

The program calculates the function 

(23) E(X) = (THS + THI )/2 - (ir/2 + Tan " 1 (dH(X)/dx) 
for many points in the interval (XSTRT, XSTOP) and when this 
function changes sign a specular point has been located. The col- 
lection of points so located is stored in an array XN(J). To save 
running time two searches are made, first a coarse grain search and 
then, in the neighborhood of each specular point, a finer grain pass 
is made. 

The search must satisfy two requirements. First, it must be 
fine enough to locate all specular points; this requires that the 
surface must be sampled often enough to get an adequate description 
of its structure. For example if the surface were described by a 
Fourier series then one would expect that sampling every twentieth 
of the minimum mechanical wavelength would be sufficient. Secondly, 
the specular positions must be located to within a small fraction of 
an electrical wavelength so that the phase relationships among the 
various specular points are correctly maintained. In the light of 


these considerations a first search might be made at a step size of 
(the minimum mechanical wavelength )/20. The fine grain search would 
then be made with a step size of say (\ e /20.0) or (1st step size/2.0) 
whichever is the smallest. In the program, the coarse step size is 
called DLTAX (DE LTA X) and the fine step size is called DLTAXOO. The 
local angle of incidence for each specular point is stored in an 
array ANG(J). This angle is used in the computation of the scat- 
tered field and is shown in Fig. 5. Once a complete pass is made 
over the surface, the scattered fields are computed. It should be 
noted that whenever any one of THI, THS, H(X) is changed, the 
complete pass must be made again. 

The actual program, given in Appendix A, makes the scattered 
field computation for two cases: 

1) all specular points contributing, 

2) scattering from concave up specular points neglected 
when calculating the scattered field. 

The second case, clearly incorrect, was an attempt to see how the 
computed fields would correspond to the results of certain statisti- 
cal theories which neglect the concave up specular points. In the 
program the electric field calculated from the first case is called 
ESCNS (ELECTRIC FIELDS SCATTERED WITH NO SHADOWING) and from the 
second case ESCDNS (ELECTRIC FIELD SCATTERED FROM CONCAVE DOWN 
POINTS WITH NO SHADOWING). 

Geometrical optics allows shadowing to be taken into account 
without much extra effort. The three types which may occur 
(specular point not illuminated by source, specular point not visible 



to observer, both) are shown in Fig. 6. Each point in the array of 

specular points, XN , is examined for inbound shadowing in the 

following way. A line is passed through the specular point XN., 

J 

H(XN.) with slope tan(THI). The equation of the line is 

J 

(24) YI(X) = Tan(THI)x + (H(XN.) - Tan(THl) XN.) 

J J 



THS 




Fig. 6. --Specular point shadowing. 



Then x is incremented in the proper direction until one of the 
following occurs. The first possibility is that at some point 
x, YI(x) becomes greater than the maximum value that H(x) can 
attain for any value of x in the interval XSTRT, XSTOP. This 
value of H(x) is called HMAX and must be supplied for each surface 
being considered. If the surface is a sum of sinusoids then HMAX 
is equal to the sum of the individual magnitudes. The second 
possibility is that at some point the value of x is incremented out 
of the interval (XSTRT, XSTOP) being considered. The third and 
final possibility is that at some point x the line YI(x) intersects 
the surface profile H(x). When the first or second case occurs the 
specular point is not shadowed. In the third case the specular 
point is inbound shadowed and for that particular j, XN( j ) is set 
equal to a number much larger than XSTOP. This allows XNj to be 
skipped when the contribution from each of the specular points is 
being computed. A very similar test is applied for outbound 
shadowing. 

When both the inbound and the outbound shadowing tests are 
completed the array of specular point positions contains values 
which are either in the range XSTRT < X < XSTOP or XNj >> XSTOP. 

The scattered field is calculated as in the case where shadowing 
is neglected except that when XN- > XSTOP the field from this 

J 

specular point is not put into the sum. The scattered field with 
shadowing accounted for is called ESCWS (ELECTRIC FIELD SCATTERED 
WITH SHADOWING) and the scattered field calculated with only con- 
cave down non-shadowed specular points contributing is called ESCD. 



C. Using the Geometrical Optics Program 

While the storage reguirement is minimal, the running time of 
this program depends largely on the step sizes which have to be 
used during the search for the specular points, and the number of 
scattering angles. This means that as the length of the surface 
increases, the time per pass required to find the specular points 
goes up and the number of passes over the surface also increases, 
since to see detail in the scattered field pattern the scattering 
angle must be examined at a larger number of points (finer grain). 
The half-power beamwidth of a uniformly illuminated aperture of 
width XSTOP-XSTRT, 

^ 0,88 x e 

(25) beamwidth = x ' STQP - ~ X STR ' T raclians 

affords a crude estimate of the fineness of the grain which must 
be taken. The increment in THS should be less than a fifth of this. 

The program has been checked for several cases, two of which 
will now be mentioned. The simplest check was the comparison with 
hand calculations for a surface described by 

(26) H(x) = 50 cos (2ttx/800) 

with x in the range (-200,200). This surface has only one specular 
point or none at all depending upon THI and THS. Another check was 
performed for a sinusoidal surface like the one shown in Fig. 7. 



y 





Fig. 7. — Specular points on a sinusoidal surface. 

In this case the specular return comes from a collection of regu 
larly spaced points which look like a pair of linear arrays of 
point sources. The program found the specular points and cal- 
culated the total scattered field correctly. 


CHAPTER III 


THE PHYSICAL OPTICS METHOD 


The next complexity of approximation to the scattered fields 
to be considered here is given by the physical optics method. 

A. The Physical Optics Approximation 

Physical optics (P.O.), (Ref. [14]), approximates the true 
surface currents on a perfectly conducting body by the currents 


(27) 


~ n 

2n x H on the portions of the surface which are 
i 1 1 umi nated 



on the portions of the surface which are 
shadowed 


where n is the outward normal to the surface and H is the incident 
magnetic field evaluated at the surface. These approximate currents 
are then used in the radiation integral to calculate the scattered 
fields. The P.0, surface current is exact if the scattering body 
is perfectly conducting half space and the incident field is a plane 
wave. As the surface curvature decreases the P.0, currents depart 
more and more from the true currents; as the curvature at some point 
on the surface goes to zero (a wedge), the method fails entirely. No 
do the scattered fields predicted by P.0, satisfy the reciprocity 
theorem except for backscatteri ng. Nevertheless, the P.0, method 
has a significant advantage over G.O. in that the fields remain 



bounded even if the radius of curvature of the surface becomes in- 
finite. Hence the flat facets of a surface can be approximately 
analyzed. 

Whether or not P.0, provides any more useful information than G.O. 
is a question of long standing, and the answer seems to depend upon the 
geometry of the scattering body (Ref. [15]). For the kind of surfaces 
considered here it will appear that P.0, gives a good approximation 
to the scattered fields over a significantly wider range of surface 
characteri sti cs than G.O. It is important to note that in this work 
the far field radiation integral over the physical optics currents is 
evaluated numerically to give the scattered fields. Unlike a number 
of rough surface scattering theories (Ref. [16]), no stationary phase 
approximation to the far field radiation integral is used. It is 
well known (Ref. [17]) that when the stationary phase approximation 
must be made, one obtains the G.O. result and there is then no dif- 
ference between the two approaches. 

The far-zone scattered fields will now be calculated using the 
physical optics currents. In the T.M. case, (see Fig. 8) the 
incident electric field is a z polarized plane wave of unit magnitude 
and the incident magnetic field is 

. +jk(xcos (THI)+H(x)sin(THI )) 

(28) H 1 = e [-si n(THI )x + cos(THI) y]/ n 

where n is the impedance of free space. Using Ref. [18] and the fact 
that the tangential electric field vanishes on the surface, the 
scattered electric field is given by 


21 




Fig. 8. --Geometry for T.M. physical optics. 


(29) 


C _ 

E ^o> = - TT 


-jk I r-r 


ifi \ e 


(nxH ) 


dz dc 


c m“°° 


r-r. 


where r Q is the position vector to the observation point, F is the 
position vector of a point on the surface and n is the unit outward 
normal to the surface. The notation c^-j indicates that the inte- 
gration is to be carried out only over those portions of the contour 

which are optically illuminated. 

— i ~ 

Since H and n are independent of z one can show, by using an 
appropriate integral representation for the Hankel function (Ref. 
[18]), that the scattered field is 


22 



(30) E s ( Pq ) = - J (2nxH^ ) H^ 2) (k|p-p o |) dc 

c in 

where all variables are confined to the x,y plane 

(31) p 0 = x 0 x + y 0 y 

(32) p" = x x + y y 

and Hj y (x) is the Hankel function of the second kind and zero order. 
Using the large argument approximation for H^(x), the far field 
scattered electric field becomes 

(33) E s z (p 0 ) = - (^) 1/2 |e J4 -? f sindHI-tan-^fi)) 

jkQ(x) j 

e Jl+(H)^ dx 

where H(x) describes the surface height profile, 

(34) A = . 

and Q(x) is given by Eq. (16). As before, the factor 
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has been suppressed in both the computed and reported values of the 

c , 

scattered electric field, so that the actual field E z (p q ) is related 
to the print out value E^ by 



When the incident magnetic field is z directed (transverse 
electric case) it is convenient to work with the scattered magnetic 
field. The latter is found from Ref. [18] 

-jk IF-7 I 

^ -i p o 

(nxH‘ ) x v — dz dc 

|r-r 0 l 

where H 1 is the incident magnetic field (see Fig. 9). The two 
dimensional far field scattering becomes from Eq. (36) 


(36) 4tt n (r Q ) = 2 



Fig. 9. --Geometry for T.E. physical optics. 



I 

jkQ(X) 


Again, the factor 

-j k |p"o I 
e _ 

Ji%i 

is suppressed in the programs of Appendix A, so that the plotted 
or tabulated field strengths, h|, are related to the true fields, 
h|(F 0 ) by 



There are two further considerations that may be discussed at 
this time. For bi static scattering it may happen that not all of 
the currents set up on the surface by the incident field are optically 
visible to the observer (see Fig. 10). In the physical optics pro- 
grams developed here no account was taken of this possibility. 

Obviously such considerations do not arise for backscattering. 

So far, in this chapter a perfectly conducting surface has 
been assumed. Physical optics can be generalized to treat dielectric 
surfaces by using a pair of equivalent electric and magnetic surface 
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Fig. 10.--0pti cally invisible surface currents. 

currents obtained from the fields of a plane wave incident on a 
dielectric half space (Ref. [19]). Since two integrations would 
be required to compute the scattered fields, it would seem that the 
running time should nearly double, but very little extra storage space 
would be required. 

B. Discussion of the Physical Optics Computer Programs 

For either polarization the physical optics program is divided 
into two parts. The first, and by far the most difficult, finds 
the shadow boundaries on the surface, since the integrations are to 
be performed only over the illuminated section of the contour. The 
second part performs the necessary integration to calculate the 
scattered far fields. 
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The program opens by considering the function H(X) which des- 
cribes the surface between the defined endpoints ALEP (Left End Point) 
and REP (Right £nd Point). The search for shadow boundaries begins 
at REP by determining whether or not the right endpoint casts a 
shadow on the surface and proceeds from right to left (see Fig. 11). 




END POINT CASTS SHADOW ONTO END POINT DOES NOT CAST A 

THE SURFACE i.e. T AN“'(j^| ) >THI SHADOW ONTO THE SURFACE 

i.e. TAN _, (^| ) < THI 
dx Lep 

Fig. 11 .--Shadowing at the right end point. 

If THI (the incidence angle-requi red to be less than 90°) is greater 
than 80° it is assumed that no shadowing occurs. The starting point 
of the illuminated zone (either REP or A of Fig. 11) is stored in the 
first position of an array called SX (Shadow boundaries X^ co- 
ordinate). The value of x is decremented until either a point on 
the surface is reached where the tangent-slope condition 

(39) tan (THI) 

is satisfied, at which point a shadow zone begins, or x becomes less 
than ALEP, in which case the second entry inSX is ALEP. On the 
other hand if Eq. (39) is satisfied for some x between SX-j and ALEP 
then this value of x is stored in SX 2 , a line with slope tan (THI) 



is passed thru the point, and its intersection (if any) with H(x) is 
found. If there are no such intersections, then all of the surface 
to the left of the point is shadowed. If an intersection does exist 
then the search for a point where the tangent-slope condition is 
satisfied begins again. This process continues until x is decre- 
mented past ALEP. The array SX thus stores the positions of 
points with an illuminated zone on their left in oddly subscripted 
locations and the points with an illuminated zone on their right in 
evenly subscripted locations (see Fig. 12). The size of the decrement 
used to locate the boundaries should be small enough to catch the 
surface features, and to locate the ends of the shadow zones within 
a fraction of a wavelength. 



Fig. 12.--I11 ustration of shadowed and illuminated zones. 

The integration over the illuminated sections of the surface 
to find the scattered fields is performed in a subroutine called 
BINT(XX,YY) (Bistatic radiation Integral) the arguments of which are 
the initial and final coordinates of one of the illuminated zones in 



sx(J). The integration is repeated for each zone until all illum- 
inated zones have been considered. The total scattered field (called 
S) for a particular THI and TNS is the sum of the zone fields. 

Except for normalization, the programs for the two polarizations 
differ only in the subroutine called FTBI(X) (Function To Be 
Integrated); the factor sin(THI-tan~^ (H)) for the T.M. polarization 
is replaced in the T.E. case by sin(THS-tan”^ (H)). The actual 
integration over the physical optics surface currents is performed by 
a five point Gaussian integration. In choosing the interval on the x 
axis over which the five point Gaussian integration is to be applied, 
two conditions must be met. The first is that the number of sample 
points along the contour must exceed five per wavelength. Presuming 
surface slopes of less than 60°, this means that ten sample points 
should be taken per electrical wavelength on the x axis. The second 
condition is that, if the surface were to be represented by a Fourier 
series, there should be 8-10 sample points per minimum mechanical 
wavelength along the x axis. Presuming s for example, that the first of 
the above conditions is the most stringent, each section of illumin- 
ated surface (i.e., between x = SX^ and x = SX^., j odd) would 
be divided into half electrical wavelength intervals plus a fractional 
interval, and five point Gaussian integration would be applied to each 
of the half electrical wavelength intervals, and to the last, 
fractional, interval . 
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C. Comnents on the Use of the Physical Optics Programs 


As in the case of 6.0. , the storage requirements are minimal, 
while running time depends upon the length of the surface and number 
of incidence and scattering angles which are investigated. For each 
THI the search for ill uni nation boundaries is performed only once, 
but the integration must be repeated for each scattering angle con- 
sidered. For many of the scattered field computations considered here 
the angle of incidence was held fixed and the scattering angle was 
varied between 0 and 180°. For such cases the time required to find 
the illuminated zones on the surface is small compared to the time 
required to do the integrations for the scattered field. 

As the surface length is increased the time required goes up 
rapidly since the integration for each scattering angle takes longer 
and THS must be incremented with a finer grain to get an accurate 
reproduction of the structure in the scattered field pattern. The 
size of the increment for THS has already been discussed in connection 
with the geometrical optics program. For example, the time required 
to run a surface 16 electrical wavelengths long, with THS incre- 
mented by 0.5° from 0 to 180°, was 1.8 min. By comparison, 21 min. 
were required for a surface 100 electrical wavelengths long with 
increments in THS of 0.25° from 30° to 170°, i.e., 560 values of THS. 
The value of the increment in the last case appears to have been 
just adequate to see the detail in the pattern. 

Among the checks of the P.0, program is a computation for a flat 
strip with no tapering of the illumination, for which a closed form 



physical optics result is easily obtained. The agreement was 
excellent for both polarizations. In Chapter V, P.0, will be compared 
with the other two methods of computing the scattered fields. Special 
attention will be given to the range of surface parameters over which 
the P.0, approximation is valid. 



CHAPTER IV 


THE INTEGRAL EQUATION METHOD 

In this chpater the third and most accurate method for calcu- 
lating the scattering will be examined. Here the scattered field 
is obtained from the exact surface current, which is found from a 
moment method solution of an integral equation (see, e.g.. Refs. [20], 
[21]). There are no restrictions on the curvature or form of the 
surface, but because of machine storage limitations only surfaces of 
rather short length (30 to 60 x ) can be handled. 

A. Moment Methods 

This section contains a brief introduction to the method of 
moments. For more information and other applications of this method 
refer to Ref. [22], on which the following is based. 

The objective of the moment method is to determine, numerically, 
the function F which is a solution of the inhomogeneous operator 
equation 

(40) C(F) = G 

where C( ) is a given linear operator and G is a given function. 
Suppose that F can be expanded in a set of basis functions b n 

N 

(41) F = Y F b 

~ n n n • 
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where is the n-th unknown coefficient of the expansion of F in 
that basis. Note that if a computer is to be used, N will have to 
be finite. Using the linearity property of C 

/ N \ N 

(42) C(F) = C [ l F b | = l F C(bJ = G. 

\n=l " "/ n=l n n 

To convert the operator equation to a set of simultaneous equations 
an inner product, a scalar, <h,g> is defined for functions h,g and s 
and scalars a, 3 such that 

(43) <h,g> = <g,h> 

(44) <ah + gg,s> = a<h,s> + 3<g,s> 

(45) <h,h*> = 0, if h = 0. 

Let {W- } be a set of weighting functions and take the inner 
product of both sides of Eq . (42) with W . Using the properties of 
the inner product, the original operator equation is converted to 

N 

(46) 7 <W , C(b )> F = <W , G> 

S m n n m 

n= 1 

which is exactly the familiar matrix equation 

(47) X fn = Gm 

where 
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( 48 ) 


C = <W m , C(b )> 
mn m n 


and 

(49) G = <W ,G>. 
mm 

The solution, F^ , to this system of equations can be found by any 
one of several methods, two of which are discussed in Appendix B. 

The solution may be exact or approximate depending upon N, , and 

v 

For the integral equations to be solved here, the current is 
expanded in a basis of non-overlapping pulses of unit amplitude, 
while the weighting functions are chosen to be delta functions whose 
singularities occur at the centers of the pulses. The inner product 
is chosen to be 


(50) 


<g,h> 


g h dc 
c 


where c is the contour of the scattering surface. This choice of 
basis and weight functions amounts to enforcing the integral equation 
at the centerpoints of the pulses, and is usually called "point- 
matching." For the operator equations considered in this work the 
system of simultaneous equations which result from point matching are 
well conditioned, i.e., suitable for computer solution (see Ref. [23]). 



B. Integral Equation for Transverse Magnetic Polarization 


In order to apply the point matching technique to the rough 
surface scattering problem, it is first necessary to find an 
appropriate linear operator. For this purpose the integral equation 
relating the unknown surface current to the known incident field 
has been chosen. 

A 

The incident electric field is z directed, the incident magnetic 
field is transverse (T.M. polarization) to the generators of the 
surface with contour c as shown in Fig. 13. If the total electric 



Fig. 13. --Geometry for T.M. scattering. 

field is written as the sum of the incident field F 1 and the scattered 
field E 5 , the boundary condition 

(51) E 1 + E 5 = 0 

must be satisfied on c. The scattered field is given in terms of 
the z directed surface currents, J z (?‘ )» by ( see Ref. [24]) 
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(52) 


E|(p)=-? n -| J z (p ' ) H^Vlp-p'D d«,‘ 

c 

(?) 

for the two dimensional case, where H' ' is the Hankel function of 

o 

the second kind and order zero, n is the impedance of free space 
and k is the wave number, 2-rr/Xg. Combining this with the boundary 
condition (Eq. (51)) gives the integral equation for the unknown 
surface current 

(53) = I 2 - | J z (^') H^ 2) (k|pV|) d*,' 

c 

where p~, p 1 are now both confined to the contour c. Equation (53) 
can now be identified with Eq. (42) as follows: 

E^(p") corresponds to G, 

J z (p"') corresponds to F, 
and the operator 

I 2 - j ( ) (k | p"~p _1 | ) d&‘ corresponds to C( ). 

c 

As it stands the integral equation requires the consideration of 
the current on the entire boundary c; if the entire contour of a 
two dimensional earth were to be included, the storage requirements 
for a moment method solution would be astronomical. It seems reason- 
able to assume that for standard radar wavelengths and with directive 
antennas, the surface current is appreciable over only a very small 
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portion of this contour. Thus it will be presumed that the surface 
current outside a certain illuminated region, which extends from -EP 
(End Point) to +EP, can be neglected (see Fig. 14). To simulate the 



Fig. 14. --Modification of true contour to 
a shortened contour. 


illumination of the surface by a directive antenna, an amplitude 
taper t(x) is introduced* in the following way. The amplitude of 
the incident field is taken as unity to within two electrical wave- 
lengths from each end point. Between one and two electrical wave- 
lengths from each end the field is sinusoidally tapered to zero. 
Over the last wavelength the incident field is taken to be zero. 

The incident field with tapering included, E^(jT), is thus 


*The use of amplitude tapering on a plane wave amounts to independently 
specifying the amplitude and phase, which makes such a field differ 
from a Maxwellian field. Changing the specification of the incident 
field would require no fundamental change in the methods and programs 
employed. 


( 54 ) 


+j|^(cos(THI) x + sin(THI) H(x)) 

4 ^ = t(x) e- jk ' p = t(x) e e 

The neglect of the surface currents beyond the endpoints (±EP) 
has been checked by lengthening the dead zone at each end of the region 
under consideration and noting the change in the surface currents and 
scattered fields. The results of this test are presented in Section 
D of the chapter and do indeed justify the assumption of negligible 
currents beyond the illuminated region. 

Although tapering of the incident field is not needed in the 
P.0, or G.O. formulations, it has usually been included in the 
calculations so that the results of all the techniques can be fairly 
compared. The only cases in which tapering is not used are special 
tests of the individual methods. 

The integral equation becomes 

EP 

(55) = Sr f J z (p "' ) H o 2 ) ( k lp'vi> de ‘ 

-EP 

with p~, p" 1 both confined to the section of the contour for which 
-EP<x<EP . 

The method of moments can now be applied. The surface is 
divided into segments of equal arclength DC, and the current, J z , is 
expanded in a basis of non-overlapping pulse functions as 
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(56) 


J z (p‘) = 


N 

l 

n=l 


F_ P. 


( p, - p n ) 


where p" n is the position vector of the midpoint of the n-th segment 
of the surface, F n is a complex number representing the magnitude 
and phase of the current over the n-th segment of the contour, and 
the n-th basis function (p"' -p" n ) is a pulse of unit amplitude and 
width DC along the contour c. Thus the actual surface current is 
to be approximated as shown in Fig. 15. For a reasonable represen- 
tation of the surface current, the pulse width, DC, must be a fraction 



Fig. 15. --Approximation of the surface current. 

of an electrical wavelength; X e /1 0 has been found to be satisfactory. 
The shape of the surface must also be considered in choosing DC, 
since the surface must be accurately modeled by strips of width DC. 
Hence, if x m is the shortest mechanical wavelength in the Fourier 
spectrum of the surface, then DC should also satisfy DC <_ A^IO. 

Of course the more restrictive of the two conditions should be met. 
Applying the method of Section A of this chapter to Eq. (55) 
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EP 

P ¥ (p'-p n ) H^ ) (k |p-p^' | ) ds.' 


Prj H o 2) ( k lp-p'|) da' 


H^ 2 ) (k|p-p'|) di' 


where / means "integrate over the n-th segment of the contour". 
DC n 

Taking the inner product of Eq. (57) with the weighting functions. 


(58) 

<«(p-P m )» E’(p)> = I 2 - 

. 1 , p n ^(P-Pm 5 • }h< 2) 

oc n 

so 



(59) 

'Xl-SM/nl 

n i nr 

N^MPm-ri) d *' 


DC n 


which is the same as the NXN matrix form 

(60) [C] [F] = [E] 

where 


(57) 


E z(p) 


kn. 

4 


-EP 


l F n 

n=l n 


- kn 


Y f 
4 L k 

n=l 


EP 

r 


-EP 


kn N f 


DC. 


( 61 ) 






d£' , 


J 
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(62) 


E = 
m 


E 1 (p ) 
z VM m' 


and F n is the unknown amplitude and phase of the current in the 
n-th contour segment. Once Eq. (60) is solved, the surface current 
is known. 

The far field scattering from the surface is found from the 
surface currents and Eq. (52) to be 




e -j' k ip| 

W\ 


EP 

■ 

-EP 


J z (p"') e jk(p,-p) dd' 


^ kri 


H P J 4 

4 Trk e 


. 5 7T .. 
3 t - 


i-i EP 




-EP 


I F n (p'-pJ e jk(p, ‘ p) dd' 
n=1 n * n 


= !SH 2 
4 \ rrk 




DC 


N 



p e jk (p *p) 


The output of the computer programs is a normalized scattered field, 
E^, which is related to the true scattered field, Eq. (63), by 


(64) E| = (p") e j k l p l . 


C . D iscussion of the Computer Program for Transverse 

Magneti c Pol ari zati on 

Several different programs were written using the above formu- 
lation of the problem. In the first part of this section the common 
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features of the programs will be discussed and later their dif- 
ferences and relative merits. 

All of the T.M.I.E. (transverse magnetic integral equation) 
programs require that the surface have its arclength subdivided into 
segments of width DC, and have the endpoints and midpoints of these 
segments stored. The surface breakdown is shown in Fig. 16. The 



Fig. 16. --Breakdown of surface into segments 
of length DC. 

3 -th segment lies between and , while the j-th midpoint 
(XM.) is such that X . <XM .<X . +1 . The surface is segmented by using 
the arclength formula and rectangular rule integration. After the 
surface subdivision is completed the programs differ somewhat 
depending on how the matrix elements are calculated. 


Once the matrix elements have been calculated the first part 
of a two part solution of the system of equations begins. In all 
of the solution methods used the matrix is factored into an upper 



and a lower triangular matrix, see Appendix B. The matrix elements 
depend only upon the surface profile H(x), and are independent 
of the incident field, THI or THS so that the factorization need be 
done only once for a given profile. In the second part of the 
solution the array [F] is loaded with the tapered incident electric 
field at each of the XM. ; the back substitutions (described in 

J 

Appendix B) are then carried out to find the current coefficients, 

F . The scattered fields are then calculated from Eqs. (63) and 
(64). 

The differences in the several programs for the T.M.I.E. lie 
mainly in the calculation of the matrix elements (Eq. (61)). The 
simplest way to evaluate Eq. (61) for n^n is to presume that 
H o 2)( k|P m -p' , l) constant over the n ~th interval; then 

< 65 > c w-M 2) < k IF n -F„l> DC 

( 2 ) 

If m=n, a small argument approximation to hP '(X) is made and 
integrated analytically, giving 

(66) C = jH- DC 

v 7 nm 4 o 2 e 7 

where e is the base of the natural logarithm. In practice the 

kn 

matrix elements are simply the Hankel function and the • DC is 
accounted for when the fields are printed out. This approximation 
results in a symnetric matrix which, if efficiently stored, requires 


only N(N+l)/2 storage locations. The length of surface which can be 

treated is increased by a factor of fE over that which can be treated 

by methods requiring the storage of the full matrix. Appendix B 

gives the details of the storage and solution methods. 

In another program, 5 point Gaussian integration. Ref. [25], is 

used to evaluate the C „ for m^n, and when m=n Eq. (66) is used. 

mn ^ 

2 

The matrix is no longer symmetric so all N terms must be stored. 

A third program was written which takes advantage of the fact 
that the currents are continuous on the surface except at sharp 
edges (Ref. [26]). Since the column vector [F] of Eq. (60) represents 
the current, continuity requires that adjacent entries be similar. 
Hence it is possible to interpolate. The currents at the even 
numbered stations may be approximated in terms of the adjacent 
currents by 

< 67 > F 2n " < F 2n-l + W 72 ' 

For simplicity, the original matrix will be assumed to be of odd 
order 

(68) N = 2 kk + 1 . 

If, for example, N=7 then, using Eq. (67) in Eq. (60), one obtains 
the reduced system 



(69) 


E, ■ C„F, ♦ ^(F,*F 3 ) ♦ C 13 F 3 + ^(F 3+ F 5 ) ♦ C 15 F 5 ♦ ^(F 5+ F 7 ) ♦ C,^ 


E 3 - C 31 F 1 + + C 33 F 3 + -f-( F 3 +F 6 ) + C 35 F 5 + -f (F 5 +F 7> + C 37 F 7 


'34 


36, 


E 5 = C 51 F 1 + “T* F 1 +I V + C 53 F 3 + "f L(F 3 +F 5 ) + C 55 F 5 + ’? i(F 5 +F 7 ) + C 57 F 7 


'54 


56, 


c c c 

E 7 = C 71 F 1 + _ ?' (F 1 +F 3 ) + C 73 F 3 + _ T^ F 3 +F 5 ) + C 75 F 5 + + C 77 F 7 


where only odd rows have been retained, i.e., Fg, F^, F^ are considered 
known. Collecting terms. 


(7°) E k - (C k] + -£■) 


F + (-*2 
F 1 1 2 


+ C k.3 + 


'k4 


)F, + {■ 


'k4 


+ C 


'k6 


k5 


)F f 


+ (■ 


'k6 


+ C k7^ F 7 


for k = 1 ,3,5,7, 

and the number of unknowns has been reduced to kk. Since matrix 
manipulations are made using regular subscripts in the machine, it is 
very desirable to relabel the coefficients in the reduced system as 
foil ows 


(71) 


o, _ C (2m-l),(2i-2) , r x L ( 2m-l),(2i) 

C mi = 2 + C (2m-1 ) , (2i -1 ) + 2 


for the “interior" columns where m=l ,2 ,3, • * • ,kk and i-2 ,3, • * • ,kk-l . 
The first and last columns of the reduced matrix are 


(72) 


C m2 " C (2m-1 ) ,1 


C (2m-1 ) ,2 

2 


m=l ,2 ,3 , • • • ,kk 
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(73) 


(2m-l),(2kk-U. 


c , _ C (2m-1 ) , (2kk-2) , c 
m,kk 2 

The C. ■ are the elements of the original NXN matrix while C! . are ele- 
1 J * J 

ments of the kkXkk reduced matrix. In the computer program the C! . are 

' v) 

called C. . while the original matrix elements C. . are labeled CO... 
i J ^ J * J 

When using the interpolation technique the surface is subdivided 
as usual except that, if an even number of segments is produced, then 
the last segment is dropped to make N odd. The system of equations 
is now 

(74) [C][FP] = [E] 

where [E] is filled with the incident electric field at the midpoints 
of the segments with odd subscripts and the matrix [C] is loaded 
according to Eqs. (71), (72) and (73). After the solution has been 
found the column vector FP ( J ) contains the currents on the segments 
with odd subscripts. The complete set of surface currents [F] is 
obtained by interpolation with 

F ? . , = FP. for j = 1 ,2 , • • • ,kk 

(75) d3 ~' 3 

F 2j = (FP j + FP j+1 )/2 for 3 = 1 >2, • • • ,kk-1 . 

Once the column vector [F] has been filled in, the calculation of the 
scattered field proceeds as in Eqs. (63) and (64). The interpolation 
technique has been applied to the program which uses Gaussian 
integration to calculate the matrix elements. 


46 



The big advantage of interpolation is the dramatic increase in 
the size of the surface which can be handled for a given storage 
capacity. If the machine can handle an arclength of L using the 
non-symmetric, non -interpolation program then the symmetric matrix 
program can handle an arclength of JIT L while the interpolation 
technique will do an arclength of 2 L with the same amount of storage. 
The interpolation program still requires that all of the original 
matrix elements be evaluated to fill in the reduced matrix (Eqs. (71), 
(72) and (73)). 

The integral equation programs require large amounts of storage 
and fairly long running times compared to either the 6.0. or P.0, 
programs. The IBM 360-75 used here can hold a 275 x 275 complex 
matrix in high speed storage so that surfaces of length 27 A., or 54 A 
if interpolation is used, can be handled with DC = A e /10. As for the 
running time, consider the 16 A g long surface mentioned in Chapter III 
Section C, v/hich took 1.8 minutes using the P.0, program. The scat- 
tering from the same surface was computed by the three T.M. integral 
equation methods. The symmetric formulation required 2.8 minutes and 
storage for 14,000 complex numbers. The program which uses Gaussian 
integration to evaluate the matrix coefficients required 5.0 minutes 
and twice as much storage^while the interpolation program required 
3.3 minutes and storage for 7,000 complex numbers. Where speed is 
important the use of the symmetric I.E. program is indicated, while 
long surfaces are best handled by the two point interpolation program. 



D. Tests of the Transverse Magnetic Integral 

Equation Programs 

The shortened contour assumption is one of the most crucial 
in the construction of the integral equation programs (Fig. 14). The 
obvious way to test it is to extend the non-illuminated portion of the 
surface, which amounts to lengthening the contour without changing 
the non -zero portion of the illumination (see Fig. 17). If the 
approximation is indeed valid, then the current in the non -illuminated 
sections should fall off rapidly and the scattered fields should be 
the same in both cases. The assumption was tested on a sinusoidal 
surface, using the program with Gaussian integration. When regular 
tapering was used, the current at the outer ends of the dead zones 
was down by a factor of 30 from that in the central part of the 
contour. When the extended taper was used, the current at the new 
outer ends was down by a factor of 100, The scattered fields for the 
two cases are displayed in Fig. 18 and show clearly that the dif- 
ferences are insignificant. Thus it may be concluded that tapering 
of the incident field does permit the replacement of the true contour 
by the shortened contour. 

The wedge. Fig. 19, for which asymptotic solutions are available, 
provides a test case for the integral equation programs. The angle of 
incidence, THI, was chosen to be 90°. In order to emphasize the 
corner contribution, a Gaussian tapering of the incident field was 
used, i .e . , 


-(x/2a )' 
t (X) = e e 


( 76 ) 



t(x) NON ILLUMINATED 

*\ /ZONE 


I 
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Fig. 17. —Contour and tapering function used to test the shortened contour assumption. 


REGULAR BOUNDARY 
EXTENDED BOUNDARY 
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THS— THE SCATTERING ANGLE 



Fig. 19. --Geometry for wedge test. 

The surface current. Fig. 20, shows the expected singularity at the 
corner. The computed scattered field is plotted in Fig. 21 along with 
the scattered field calculated independently using the geometrical 
theory of diffraction, Ref. [27], Again, the agreement is seen to 
be excellent. All three T.M. integral equation programs produced 
essentially identical scattered fields. In a test of the self 
consistency of the three programs the scattering from the surface 

o 

H ( X ) = 5 sin x was computed. The differences in the scattered 
fields are very minor and would not be perceptible on the scale of, 
e.g.. Fig. 18. 

In the light of the above tests, there seems to be no reason 
to prefer one T.M. integral equation program over the other two if 
numerical accuracy is the only criterion. If the running time or 
storage requirements must be considered then the preferred formulation 
can be determined by the comments at the end of Section C of this 
Chapter. 
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SIMPLE DIFFRACTION 


Fig. 20 .--Computed |J 



on a wedge, T.M. case 
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E. Integral Equation for Transverse 
Electric Polarization 

For the T.E. polarization, the incident magnetic field FT 1 is 
z directed and it will be convenient to work with the integral 
equation for the magnetic field given (Ref. [28]) by 

-jk|F-r‘ | 

(77) T (r) = 2 n x FT 1 (r) + n(r) x f J _(r' ) x V 1 ds 

S C.TI J S |— _ r , | 

where F, F 1 are both position vectors of points on the surface, FT 1 (r) 

_____ A 

is the incident magnetic field, J s (r) is the surface current, n is the 

outward normal to the surface and j indicates that the region about 
— — s 

r‘ = r is to be deleted from the integration. See Fig. 22. 


A 



Fig. 22. —Three dimensional geometry for 
T.E. integral equation. 

The two dimensional integral equation can be obtained by con- 


sidering an infinitely long cylinder as shown in Fig. 23. When the 
incidence direction lies in the x,y plane the fields and surface 
current have no z dependence so that Eq. (77) can be reduced to 
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x,y 


PLANE 


Fig. 23.— Two dimensional geometry for 
T.E. integral equation. 

(78) T $ (p ) = 2 n^T) x H* Q + |jn(p) x | J^') x (r? ) 

c 

H{ 2) (k|F-7D dc' 

where (p-p ' ) is the unit vector in the p-p* direction and Hj ; (x) 
is the Hankel function of the second kind and order 1. 

Just as in the T.M. case, tapering is introduced to account 
for the directional properties of radar antennas, and to limit the 
size of the system of linear equations which will result from 
Eq. (78). One may now assume that the surface currents are zero 
except near the illuminated region and the closed contour can then 
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be replaced by the open contour of Fig. 24. For this polarization 
the current flows transverse to z along the surface so 



Fig. 24. — Open contour. 

(79) J s (p') = U x n (^)) J s (p') = T(i?) 0 s (p') 

A __ ✓V _ 

where T(p') and n(p') are the unit tangent vector and the unit normal 
vector to the surface, as shown in Fig. 24. T(p') is given in terms 
of the profile, H(x), by 



where H has the meaning assigned by Eq. (34). Using 
(81) dc' = (1 + (H(x' )) 2 ) 1/2 dx' 
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and Eqs. (78) and (79) with the tapered incident field 


(82) 


-jk, *p 

(p) = t (x) e 


the integral equation becomes 


(83) 


-t (x) 


jk-j -P 


j s (p) 

2 


+ 


jk 

4 



J s (p')H 1 (2) (k|p-p , |) 

v (x-x') 2 +(H(x)-H(x')) 2 


‘[(H(x)-H(x' ))-H(x)(x-x’ )] dx‘ 

where the integration over x' excludes a small region in the contour 
about the point described by p*. 

The mathod of moments is applied to Eqs. (83) just as in the 
T.M. case, the current is expanded in a basis of non -overlapping 
pulse functions of width DC, delta functions are used as weighting 
functions and the scalar product is the same as in the T.M. case. 

The current is thus represented by 

N 

(84) J s ( P ') = l F n (p'-p n ) 

where, p"' , p n lie on the contour c and p n is the position vector of 
the midpoint of the n-th segment, the F n 's are the unknown expansion 
coefficients and the pulse functions (pT'-fT ) have been described 
in connection with the T.M. case. Placing this current in Eq. (83), 
taking the scalar product of both sides with the weighting functions 
and using the non-overlapping property of the basis functions results 
in 
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I 


(85) 


-jk,- Pm F 

-t(x m ) e 1 ro = / + 


nr 


EP 


^ I F n | P ¥ * k l p m‘ p '^ 

n=1 ip 


( 2 ) 


[(H(x m )-H(x , ))-ft(x m )(x m -x 1 )] 
v (x m -x , ) 2 + ( H(x m )-H(x')) 2 


Since it is necessary to avoid p‘ = p m in the integration of Eq. (85), 
the summation will be forced to skip n=m giving as a system of 
equations 


( 86 ) 


-jk.*p m N 

-t ( x ) e = l C mn F 

' m mn n 

n=l 


nT 


where 


(87) C 

mn 1 


2 if rn=n 


*n+l 


, k r / p \ L(H(x )-H(x'))-H(x )(x -x’ ) 3 

4 1 H { (k l p m -P' 1 > ===== — I*' 


J 

v 'Si 


,(x ra -x') 2 + (H(x m )-H(x')) 2 


i f m^n 


and x n+ -| , x n are the upper and lower x coordinates of the endpoints 
of the n-th surface segment respectively. 
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Once Eq, (86) is solved for the coefficients of the surface 


current, F , the scattered field may be found from Ref. [18] 
m 

-jk|r-P 1 

(88) ^(F) = f J (?') xv' ^ ds '. 

4,r J s | r-r 1 1 

Specializing this to the far field scattering from an infinite 
cylinder and using the fact that T (p ) is independent of z and 
non zero only over a portion of the cylinder (see Fig. 25). 


/ 

/ 



Fig. 25. --Geometry for calculation of far 
field scattering, T.E. case. 

-jk|p] i^L EP 

( 89 ) H s (? » § i e 4 [ J (?) L s -in (THS)-H(x' ) cos_(THS] 

J]?" -EP S Jl+(H(x')) 2 

e jk(x'cos(IHS)+H(x' )sin (THS) ) dc , 


Substituting Eq. (84) whose coefficients are now known into Eq. (89) 
and assuming that the integrand is nearly constant over a surface 
segment of length DC, 


(90) 




2vT 


DC 


-jk | p | 

e 



N 

l F n cos(THS-THN(XM n )) 


jk ( XM n cos (THS )+H (XM n )sin (THS ) ) 
• e 


where THN(x) (JHETA NORMAL) is given by 


(91) THN(x) = U/2) + tan" 1 (H(x)) 


as shown in Fig. 25. The computed and plotted value of the scattered 
field, H^, is given by 

c c _ I + Jk | p~| 

(92) h| = H|( P ) Jlfle 

F. Discussion of the Computer Program for the 
Transverse Electric Polarization 

The programs for the T.E. polarization are very similar to those 
for the T.M. polarization. As in the T.M. case the contour is broken 
up into segments of equal length DC. The same notation is used for 
the endpoints (x) and midpoints (XM) of the segments (Fig. 16). The 
T.E. and T.M. programs differ mainly in the values of the elements 
of the matrix [C], and in the driving side of the system of equations. 


Also, for the integral equation used, the matrix is n on -sy mne trie no 
matter how the coefficients are evaluated. Once again the system of 
equations, (Eq. (86)), is solved in such a way that different 
scattering and incidence angles do not require a completely new 
solution. Only the back substitution portion need be repeated (see 
Appendix B). 

Several different programs have been written for the T.E. case, 
the major difference between them being the method used to evaluate 
the coefficients (Eq. (87)). The simplest way is to assume that 
the integrand is constant over the strip width so that 


(93) 


C 

mn 1 


2 if m=n 


jk 


(DC) 


H l 2) < k |p m -P n l> 


[(H(XM m )-H(XM n ))-H(XM m ) 


p m" p n 


(XM -XM )] if n itfn. 
m n' J 


In practice, only the five point Gaussian integration was used to 
evaluate the off diagonal elements of [C], since it did not require 
much more running time than the simpler method. However, the inter- 
polation technique retains all of its advantages and goes exactly as 

in the T.M. case with the C! . given by Eqs. (71), (72), and (73). 

* J 

Thus surface lengths of 27 a_ (or 54 A_ with interpolation) can be 
handled. As an example of the running times required, consider again 
the surface of length 1 6 A e mentioned in Chapter 3 Section C. The 
T.E. physical optics program required 1.8 minutes while an equivalent 
run using the T.E. integral equation program required 5.0 minutes. 



The interpolation program for this polarization took 3.5 minutes. 

Thus the interpolation program is superior to the non -interpolation 
program both with respect to storage requirement and running time. 

G. Tests of the Transverse Electric Integral 

Equation Programes 

The shortened contour assumption plays the same role and is 
tested in the same way in the T.E. integral equation programs as 
in the T.M. case. The contour is extended as shown in Fig. 17. 

When the regular tapering was used, the current at the outer ends of 
the dead zones was down by a factor of 70 from that in the central 
portion of the contour. When the extended surface was considered the 
current at the new outer ends was down by slightly more. The nearly 
identical scattered fields for the two cases are shown in Fig. 26. 

The wedge provides a test case for which an independent result 
is available. The test geometry is as shown in Fig. 19 except that 
here the incident magnetic field is parallel to the comer of the 
wedge. Gaussian tapering of the incident field, Eq. (76), is used. 

In contrast to the current singularity in the T.M. case, the surface 
current in the T.E. case. Fig. 27, shows the expected r^^ behavior 
at the corner. The excellent agreement between the scattered fields 
calculated by the integral equation method and the fields obtained 
from the geometrical theory of diffraction, Ref. [ 27] » is illustrated 
in Fig. 28. Both the non-interpolation and the interpolation T.E. 
integral equation programs gave the same result in this test. 
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REGULAR BOUNDARY 
EXTENDED BOUNDARY 




SIMPLE DIFFRACTION 
COEFFICIENT FAILS 


Fig. 27.--Conputed |J | near corner of wedge, T.E. case 


INTEGRAL EQUATION 




SIMPLE DIFFRACTION 
COEFFICIENT FAILS 




The consistency of the two T.E. integral equation programs was 
checked on a surface with a height profile H(x) = 5 sin(2Trx/200). The 
results were nearly identical. 

The above tests indicate that so far as numerical accuracy 
is concerned the non-interpolation and interpolation T.E. integral 
equation programs do not differ. The interpolation program is pre- 
ferred however because of the savings in storage. 



CHAPTER V 


APPLICATIONS 


In this chapter the previously developed computer programs will 
be used to check the applicability of the geometrical optics, physical 
optics and perturbation approximations to the calculation of the 
scattering from non-uniform surfaces. The integral equation programs, 
which are believed to be exact, are used as standards. 

The first surface to be considered has been especially chosen 
so that it fulfills the requirements necessary in order that physical 
and geometrical optics both give a valid approximation to the true 
scattered fields. The surface, a single half-cycle of a sine wave, 
has a profile H(X) = 50 cos ( 2 ttX/ 800 ) with x between 200.0 cm and 
-200.0 cm, and clearly has but one specular point. The incident field 
is tapered, and has an electrical wavelength of 25 cm. Unless other- 
wise noted, these conventions have been used throughout. The criteria 
for the successful application of G.0. and P.0, are met by this profile 
since the minimum radius of curvature is 12.8 and, having a maximum 
height of two x , there are several Fresnel zones on the surface. The 
scattered fields predicted by the G.O., P.0, and I.E. programs are shown 
in Figs. 29 and 30 for the T.M. and T.E. polarizations respectively. 
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It is apparent that all methods give nearly the same result for 
THS between 87° and 155°. No scattered fields are predicted by 6.0. 
for THS outside the range 78° and 163° since the normals to the surface 
have a limited range of directions as illustrated in Fig. 31. The 



Fig. 31 .--Limi tation of scattering directions predicted 
by geometrical optics. 

rise in the value of scattered field predicted by G.O. near 78° and 
163° is due to the movement of the specular point into a region of 
the surface of increasing radius of curvature. However, as the 
specular point gets within two wavelengths of either endpoint the 
tapering of the incident field suppresses the expected singularity 
in the scattered field. 

It should also be noted that for the P.0, results, the T.M. 
fields differ slightly from the correct fields for THS near grazing. 



For either polarization the ripple observed in the scattered field 
and correctly predicted by P.0, is probably a consequence of the 
finite length of the surface. G.O., being a purely local theory, 
will not predict effects of this nature. 

As a further check of the programs, the above profile was 
multiplied by minus one, i.e., instead of being concave down the 
surface was concave up. The amplitudes of the scattered fields re- 
mained unchanged but they all showed a phase shift of 90° due to what 
in G.O. theory is termed the caustic correction factor. 

In order to establish more quantitatively the limitations on 
the G.O. and P.0, approximations, the scattered fields have been 
computed for a set of surfaces with height profile 

(94) H(X) = A sin(27rX/200) -200 cm. < x < 200 cm., 

i.e., the surfaces are two complete mechanical wavelengths long. 

With THI fixed at 60°, the amplitude. A, was varied over a range of 

5.0 cm. to 50.0 cm. so that the minimum radius of curvature, r , 

cm 

varied from 8.0 A g to 0.8 x . The important features of the scattered 
fields over this range of r cm for each polarization are shown in 
Figs. 32-37 in order of decreasing r Qm . Some general trends are 
worthy of mantion. 

In the first place, as r cm A e decreases from 8 to 0.8, the 
agreement between the P.0, results and the exact fields goes from 
excellent to poor. It would appear that as long as the surface al- 
ways has r /x_ greater than, say, 2.5, the P.0, approximation will 
cm s 
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Fig. 35. — 'Scattered fields predicted by P.O., 6.0. , and I.E. 

methods for H(x)= 15 sin(2nx/200) , T.E. polarization 





■ 60 * 



THS-THE SCATTERING ANGLE 

Fig. 36. — Scattered fields predicted by P.O., G.O., and I.E. 

methods for H(x)= 25 sin(2*x/200) , T.M. polarization 


MINIMUM RADIUS OF CURVATURE = r cm = 1 .66 X 



THS-THE SCATTERING ANGLE 

Fig. 37. — Scattered fields predicted by P.O., G.O., and I.E. 

methods for H(x)= 25 s i n ( 2 tt x/ 200 ) , T.E. polarization 


give reliable values for the scattered field. Even for values of 
r /X = 1, P.0, may still be considered usable, that is, it will 
reproduce the general structure of the scattered fields although with 
significantly lower accuracy. This limitation on the radius of curva- 
ture necessary for the successful application of the P.0, approxi- 
mation is in agreement with the results of Ref. [ 29l in which the 
current on a sinusoidal surface of infinite extent is found. Except 
for scattering and incidence angles for which no specular points occur 
or for which a specular point coincides with a point of infinite 
radius of curvature, the 6.0. and P.0, approximations give scattered 
fields very similar to each other even when they are not correct, e.g.. 
Fig. 38. It is interesting to note that where the I.E. and P.0, (and 
hence the 6.0.) fields agree the T.E. and T.M. fields are nearly 
identical but as the radius of curvature decreases the exact fields, 
T.E. and T.M., not only differ from the respective P.0, fields but from 
each other. This behavior is not entirely unexpected since for bodies 
with large radius of curvature in terms of wavelength the polarization 
independent 6.0. is known to be a good approximation. As the radius 
of curvature goes to zero, e.g. a wedge, 6.0. and P.0, both fail and 
the scattering is polarization dependent (see the wedge tests in 
Chapter IV). 

The failure of 6.0. when no specular point occurs on the surface 
or when a specular point coincides with a point of infinite radius of 
curvature makes it far less attractive than P.O., especially when 
numerical methods are involved. For example, when A=5, (see Fig. 32) 
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G.O. predicts no scattered field outside the range 102° < THS < 138°, 
and gives fields which are singular at either end of the range. On 
the other hand, the P.0, approximation correctly predicts the scattered 
fields for a far wider range of THS, including backscatter, and the 
fields are always bounded. 

It is also of interest to note that what might be called the 
"fine structure" of the scattering, particularly for THS < 80°, 

(see Fig. 32) is not due entirely to the finite length of the 
illuminated region as in Figs. 29 and 30 but is strongly controlled 
by the height profile. 

Another approximate theory whose validity can be checked by 
the numerical mathods developed here is the perturbation theory for 
the scattering from "slightly rough" surfaces as formulated in 
Refs. [ 30] and [31 ]. Perturbation theory predicts that if the ampli- 
tude of the surface profile is much less than the electrical wave- 
length of the incident fields, then the amplitude of the scattered 
field due to the perturbation of the surface is proportional to the 
surface height amplitude. This was checked by calculating, using 
the T.M. integral equation program, the scattering from a surface 
profile described by 

(95) H(x) = c (sin(2Trx/50) + 1/2 sin(2-rrx/19.71 ) ) 

for various values of c. The field scattered by slightly rough 
surfaces is dominated by the scattered field from the unperturbed 
surface (c=0) which is quite complex for the finite strips considered 



here. Thus the behavior of the perturbed fields can best be 
illustrated by considering the difference between the actual field 
and the flat plate field. The perturbation in the scattered field, 
Ep, due to the perturbation in the height profile of the originally 
flat strip is then given by 


(96) 


E = E? - E s 


zo 


where is the total scattered field as predicted by the computer 
program, and E^ 0 is the field scattered when c is zero (i.e., a flat 
strip). In order to test the prediction that |E | ac, a low value 
of c (c=0.01 cm.), was chosen as a reference surface amplitude with 
reference scattered field |E -j|, so that for a fixed scattering angle 


(97) 


'PI 1 


expresses the perturbation theory result. The exact fields are 
compared with perturbation theory in Fig. 39 for several values of c. 
The theory appears to fail at about c/c-j = 200 which corresponds to 
a root mean square surface amplitude of approximately x / 1 0 . 

In addition to permitting the examination of the applicability 
of various electromagnetic approximations to the ocean surface scat- 
tering problem, the programs permit direct calculation of the scattered 
fields from any appropriate surface. One such application is to the 
calculation of the expected value of the backscattered power from an 
ensemble of ocean-like surfaces. Such an ensemble may be constructed 
from the known height spectrum. Ref. [ 32 ]. For the sea surface, the 



c, ( c, sin ( 2-irx/ 50) + Y sin ( 2irx/l9.7l )) 
TH I =60° 

THS = 90° 

FIELD RATIO PREDICTED 

BY INTEGRAL EQUATION 

FIELD RATIO PREDICTED 

BY PERTURBATION THEORY / 
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height spectrum. Fig. 40, decays approximately as k m over the signifi- 
cant range of wave numbers, where k^ is the mechanical wavenumber. 



Fig. 40. --Sea surface height spectrum. 

Thus a particular member of the ensemble can be chosen to be a finite 

sum of sinusoids with random phases whose amplitudes vary roughly as 
_2 

k . If the k 's are not harmonically related, the surface, like the 
ocean, will be aperiodic. One example of a surface of this type is 
given by the series 

(98) H(x) = 2. 5(0. 4 sin (2ttx/200 .0 + 0.78) 

+ 0.8(10. 0/20. 0) 2 si n(2nx/l 0.954 + 1.6) 

+ 0.8(6. 66/20. 0) 2 sin(2irx/6.28318 + 2.4) 

+ 0.8(5. 0/20. O) 2 sin(2nx/4.795 + 0.4)) 

illustrated in Fig. 41. An ensemble of surfaces of finite length can 
be generated by using successive non-overlapping sections of this 
surface . 
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Fig. 41. --Four component representation 
of the surface 

Physical optics was used to calculate the expected value of 
the backscattered power and field strength from a 75 member ensemble 
made from the surface described by Eq. (98). Each member of the 
ensemble was 75 electrical wavelengths long. On a CDC 6600 computer, 


i 
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the time required for the run was about 40 minutes. The expected 

s 2 s 

values <|E|| > are shown in Fig. 42; the expected value of E z was 
found to be extremely small compared to the root mean square field. 



Fig. 42. --Expected value of backscattered 
| E| |2 from ensemble. 

Notice that no special form of the slope distribution or 
other statistical properties of the surface have to be assumed. It 
is also possible to use a point by point, i.e. discrete, representa- 
tion of the surface, such as might be generated by the prescribed 
statistical properties of the surface. 


CHAPTER VI 


SUMMARY AND CONCLUSIONS 


In this work the scattering properties of cylindrical 
rough surfaces have been investigated by several numerical techniques 
in order to test the validity of previous theoretical work. The 
results, using as checks the integral equation solutions, show that 
geometrical optics is not usable for surfaces with radius of 
curvature smaller than 2.5 a 0 and may give poor results even when 
this condition is satisfied should the scattering geometry be such 
that no specular point exists or a specular point coincides with a 
point of infinite radius of curvature. With the exception of these 
two cases, geometrical optics and physical optics give nearly the 
same scattered fields. 

It was found that the numerical evaluation of the scattered 
fields from the physical optics currents gives good results for 
almost any geometry (except perhaps deeply shadowed configurations) 
as long as the radius of curvature condition, r > 2.5 a , is 
satisfied. Physical optics, although not always so accurate, has 
an advantage over the integral equation formulation in that the length 
of surface which can be treated is not limited by machine storage 
capacity. 
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The integral equation program has been used to check the pert- 
urbation theory prediction that the amplitude of the scattered field 
increases in proportion to the increase in the amplitude of the 
surface height profile. The numerical results confirm in a quantita- 
tive way the fact that the theory fails when the root mean square 
height is about one tenth of an electrical wavelength. 

The physical optics program, because of its ability to handle 
long surfaces and its superiority to geometrical optics, has been 
applied to the direct calculation of the expected value of the 
scattered power from an ensemble of ocean-like surfaces which were 
constructed from a height spectrum similar to that of the sea. The 
computer time required, while lengthy, was not found to be prohibitive. 

The extension of the programs to very long surfaces, to non- 
cylindrical surfaces or to dielectric surfaces appears feasible only 
for the G.O. and P.0, methods; the storage requirements for an I.E. 
solution in either case would be prohibitive. P.0, would probably be 
the easiest to modify to non-cylindrical surfaces, especially if 
shadowing were neglected. Since location of the specular points becomes 
much more complicated in the non-cylindrical case, the G.O. method 
would be more difficult to implement. 
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APPENDIX A 


COMPUTER PROGRAMS 

A listing of all the programs discussed in the text is presented 
here. To facilitate understanding of the programs, the symbols used 
in the programs have been used in the text whenever possible. 

All programs require the plot subroutine listed at the end. 

The function subprograms AHAN20(x) and AHAN21(x) are required in 
the T.M. and T.E. integral equation programs respecti vely . 



no n o non oooonnoooooo 


THIS PROGRAM IS FOR BISTATIC BACK SC A TTER ING 

E SONS IS THE RETURNED E FIELD WITH SHADOWING NOT ACCOUNTED FOR 
ESCWS IS E SCATTERED WITH SHADOWING ACCOUNTED FOR 

GEOMETRICAL OPTICS FOR THE OCEAN SURFACE 

SPECULAR POINT SEARCH IS DONE IN TWO STEPS 
Hi IS MECHANICAL WAVELENGTH DEPENDENT , HZ ISREFINNED MECHANICAL OR 
ELECTRICAL WHICHEVER IS MORE STRINGENT. 

ULTAX IS THE SEARCH S IZ Efcl ,DLTAXOO IS SEARCH $IZE#2 
DELSHA IS SHADUW TEST STEP SIZE 

THIS PROGRAM CAN HANDLE 2O0 SPECULAR POINTS /PASS IE. ONE THIG1THS 
DIMENSION XN (200) .ANGLE! 2(50 ) 

DIMENS ION AC DNS ( 720) , AWS( 720) ,AWCS(720) ,ASNS( 720) , A OS (7 20 ) 
DIMENSION ECDNS ( 72 0 ), E WS (720 ) , EW CS ( 7 20 ) , V( 10 ) ,ESNS ( 72 0) 

REAL PI , P I 2 
REAL MT WO 

COMPLEX ESCNS, ESCWS, ENS 

COMMON CA,CB,CKA ,CKB , PHA , PHB , CC ,C KC , PHC 
COMPLEX E SCONS, ESCD 

NAM EL I ST/ CAT/ CA .CBtCKA.CKBfPHAjPHB.CCTCKC.PHC.WAVE.THID 
NAME LI ST/CUT/ESNS, ASNS .ECDNS , ACDNS, EWS , AWS , EWCS , AWC S , AOS 

THE FUNCTION WHICH DESCRIBES THE SURFACE IS 

H(X)=CA*S1N((CKA*X)+PHA)+CB*S1N( ( CKB* X) +PH0 )+CC*SIN( (CKC*X) +PHC) 

C A= 10 .0 

CKA=6. 28318/200.0 
PHA=C .0 
CB=0 . 0 
CKB=0 • 0 
PHR=0 .0 
CC=0 .0 
C KC = 0 . 0 
PHC=0 .0 

HMAX=ABS(CA)+ABS(CB)+ABS( CC) 

P I = 3. 1 A 1 59 
P 12= l .5707963 
TP I =6.2831 85 

WMMIN IS THE MINIMUM MECHANICAL WAVELENGTH 
WMMIN=TPI/AMAX1 ( C KA , C K8 , CKC ) 

DLX=0. 01000 
TWDLX=20.0*DLX 

NANI IS THE NUMBER OF ANGLES TO BE INVESTIGATED 
NANI =360 
X S T RT =-200 . 0 
X S TCP=-XS T RT 

THS IS THE ANGLE BETWEEN THE POS. X AXIS AND THE SCATTERING DlKEC. 
THI IS THE ANGLE BETWEEN THE POS. X AXIS AND THEINC. DIRECTION 

TH I =60. 0*3. 141592 7/ 180. 0 
C WAVE IS THE ELECTRICAL WAVELENGTH 

WAV E= 2 5 • 0 
DLTAX=WMMIN/10.0 

DLTXOO=AMI NLUDLTAX/ 5.0) ,( WAVE/20.0 ) ) 

delsha^mmin/io.o 

XSKIP=XSTOp+( 1 0**9 J 
TANTHI=TAN( THI ) 

THI D=THI*1 80. 0/3. 14159 
CSThI =COS ( TH I ) 

SNTHI=SIN( THI ) 

NAMEL IST/TOM/DLTAX,DLTXOO,DELSHA 
WRI TE( 6, TOM) 

DO 93 IRE=1,NANI 
ASNS ( IRE ) =0.0 
AC0NS( IRE )=0.0 
AWS ( I RE ) =0.0 
AWCS( IRE) =0.0 
ESNS( IRE) =0. 0 
ECDNS ( IRE ) = 0 . 0 
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EMS ( I R0 I =0 • 0 
EWCSf I«E)=0.0 
93 CONTINUE 

DO 17 I J = 1 1 NANI 

THS=FLOAT( I J)*0. 8726646 E-02 

THSD=THS*57. 29578 

AOS ( 1 J ) =THSQ 

WRITEI6, 356) THID,THSD 

356 FORMAT ( 1 IH INC ANGL E= , El 5 . 8 , 1 3H SCATT ANGLE= » El 5. 8 J 
SUCOS=CSTHI+COS( THS ) 

SUSIN=SNTHI+SIN ( THS ) 

N=0 

FIRST FIND POSITICNNS CF SPECULAR RETURN AND STORE THEM 
THE FIRST POSITION CAN NOT BE A SPECULAR POINT 
XP=XSTRT 

SUMC)2=< THI + THSI/2.0 
E=SUflD2-( TH ( XP ) +P I 2 ) 

102 X P= X P + DL T AX 
E0 = E 

E=SUMD2-(TH(XP)+PI2) 

I F ( E.EQ.O«0 ) GO TO 100 

IF( ( ( EO.GT.0.0) . AND. ( E. LT.O.OI ) .OR. ( ( EO . LT.0.0 ) .AND. (E.GT .0.0 ) ) ) 
2 GO TO 100 
GO TO lOl 

100 N=N+ 1 
XN(N) =XP 

ANGLE I N) = THS-( TH ( XP) +P 12 ) 

101 IF (XP.LE. XSTOP) GO TO 102 
IF(N.EQ.O) GO TO 372 

C THIS IS TO REFINE THE POSITION OF THE SPECULAR POINT 

DO 25 K = 1 ,N 
X S0 = XN { K ) -DLTAX 
E = SUMD2-(TH(XS0>+PI2 ) 

222 XSO=XSO+DLTXCO 
E0= E 

E = SUMD 2- ( TH(XSO) +PI2 ) 

IF(E.EQ.O.O) GO TO 252 

IF ( ( { EC.GT .0.0 ) .AND. < E.LT .0.0 ) ) .OR. U EO.LT.O. 0) .AND. < E.GT. 0.0 >)) 
2 GO TO 2 52 
GO TO 253 

252 XN { K ) =XS0 

ANGLE ( K)=THS-(TH(XSOI +PI2 ) 

253 CONTINUE 

IF (XSC.LT.XN(K) ) GO TO 222 
25 CONTINUE 

ESCNS-=CMPLXCO. 0,0.0) 

ESCDNS=CMPLX(0 .0,0.0 ) 

DO 10 K=1,N 

P HAS E = ( TP I/WAVE ) * ( { SUCOS*XN( K) ) + ( SUS I N*H l XN ( K ) ) ) ) 

KC=R5(XN( K) )*COS ( ANGLE ( K ) > 

IFIRC.LT. 0.01 PHASE=PHASE+(PI/2.0) 

ENS = ~( (SQRTl ADS<RC/2.C) )) *CEXP ( CMPL X< U. O f PHASE I M 
C TAPPERING INCLUDEO 

XG= XN ( K ) 

I F ( XG.GT. (XS TUP -WAVE) ) ENS=C MPL X ( G . 0 , C . 0 ) 

I F( XG.LT. ( XSTRT+ WAVE)) ENS =CMPLX (0 .0,2 .0) 

I F t (XG.GT .(XSTGP-I2.C *WAVE) ) ) . AND . I XG. L E . ( X STOP- WA VE ))) 

2ENS=ENS* (0. 5-(0. 5*S IN < ( 3 . 14 15 9/ WAV E ) * ( X G- ( X STOP-( 1. 5* WAVE )))))> 
IF { ( XS.GE. ( X S TR T + WAVE ) ) . A ND . < XG . L E . ( X « TRT + ( 2.0*WAVE ) ) ) ) 

2ENS = ENS* (0. 5+ (0. 5*S IN ( ( 3. 14 159/ WAVE ) + t XG-< XSTRT+ ( 1. 5* WAVE I) I I ) I 
ESCNS-ESCNS+ENS 
1F( RSI XN< K) ) .LE.O.O) GO TO 10 
ESCDNS=ESCCNS+ENS 
10 CONTINUE 

ACD- CA6S( ESCDNS) 

IFIACD.LT.l .0 E-05) GO TO 59 

ANACC=5 7.29 5 78*ATAN2 ( A I MAG ( ESCDNS ) ,REAL( ESCDNS ) ) 

59 CONTINUE 

1HACC.LT. 1.0 E-05) ANACD=0.0 
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ESMAG=CABS(ESCNS) 

ESANG=ATAN2< AIMAG(ESCNS) , REAL IESCNS) ) * 1 80 .0 /3 . 14 1 592 7 
WRITE{6,726) ESMAG, ESANG 

726 FORMAK* ','MAG. OF SCATT. E FI EL D=* , El 5 .8 , * PHASOR ANGL E= • ♦ E 15 . 8, 
23X » 1 W I THGUT SHAOPWING' ) 

WRITE ( 6 f 1 2 1 ) ACOfANACC 

121 FORMAT!' • ♦'SCATT. FIELD NO SHADOW CONCAVE DOWN TIPS ONI Y= * , E 1 5. 8 , 
2 ' PHASOR ANGLE=' ,E 15. 8) 

ESNS.t I J ) = ESMAG 
ECDNS! I J) =ACD 
ASN $ ( I J ) =ESANG 
AC DNS ( I J I = AN ACC 

C NOW FIND THE SHADOWING EFFECT 

C INBOUND SHADOWING 

IF (ABS (THI-P 12 ) .LT.0.C5 ) GO TO 500 
DO 32 7 K= 1 » N 

BI=H(XN( K ) )-( TANTHI *XN(K) ) 

STEP 1 =DELSHA 

IF ( TANTHI. LT. 0.0) STE PI =-DELSHA 
X I = XN { K ) +ST E P I 
GO TO 471 
47C XI = XI + S7EPI 
471 YI = ( TAMHI*X I )+BI 

IF (YI.IE.H(XI) ) XN< K ) = XSK I P 
I F C ABS( XN! K) ) .GT.XSTOP) GO TO 499 
IF(ABS(XI ).GT.XSTOP ) GO TO 499 
IF(YI.LE.HMAX) GO TO 470 

499 CONTINUE 
327 CONTINUE 

500 CONTINUE 

C GUT HOUNO SHADOWING 

I F ( ABS ( THS-P 12 ) -LT .0 .05 ) GO TC) 639 
TAN THS= TAN( THS J 
DO. 633 KK= 1 ♦ N 

IF (XM(KK). GT.XSTOP) GO TO 633 

THE ABOVE CARD MAKES SURE ThAT TIME IS NOT SPENT ON A PT, ALREADY 
KNOWN TO BE SHAOCWED 
BO = H(XN(KK> )-(TAMTHS*XN(KK ) ) 

STEPO^OELSHA 

IF(TANTHS.L T.O .0 ) S T EPO=-DELSHA 
XO=XN< KK) +STEPO 
GO TO 671 

670 XO=XO+STEPO 

671 Y0=( TANTHS^XQ) +B0 

IF ( YO.LE.H! XO) I XN(KK)=XSKIP 
IF(ABS(XN(KK )). GT.XSTOP) GO TO 699 
IF { A BS ( XC) ) . GE . X S TOP ) GO TO 699 
IF (YO.LE.HMAX) GO TO 670 
699 CONTINUE 
633 CONTINUE 
639 CONTINUE 

C END OF SHADOWING EFFECT 

I Nl Nl N=0 

ESCWS = CMPLX(0. 0,0.0) 



ESCD=CMPLX( 0.0*0. 0) 

DU 19 K= 1 , N 

C NEXT CARD SKIPS THE SHADOWED SPECULAR POINTS 

IF (XN(K > .GT.XSTflP ) GO TO 19 
I NI N I N=K 

PHASE=< TPl/ WAVE)* ! I SUCOS*XM K) ) + ( SUS I N*H IXN ( K > ) ) ) 

RC=RSI XN(K) )*COS( ANGLE (K) ) 

I FtRC.LT.O.O ) PHASE=PHASE+(PI/2.0 > 

ENS=-( ISQRTI A8SIRC/2.0) ) ) *C £XP{ CMPLX 1 0 . 0 , PH AS E) ) ) 

C TAPPERING INCLUDED 

XG=XN( K) 

I F( XG. GT . I XS TOP-WAVE) ) EN S=C MPLX (0.0 ,0.0 ) 

IFIXG.LT. (XSTRT+WAVE) ) ENS=C MPL X I C .L , 0 . 0 ) 

I F( ( XG.GT .{ X STOP- (2.3* WAVE) ) ) . AND .( X G. L E .( XS TOP-WAV E ) ) ) 
2ENS=ENS*(G.5-(0.5*S INI I 3. 14159/WAVE >*{ XG-I XSTOP-I 1 . 5* WA VE ) ) ) ) ) ) 

IF! ( XG.GE.I XSTRT+WAVE )).AND.(XG.LE.( X ST RT + ( 2 . C *W A VE > ) ) ) 

2EN$=ENS* I C. 5+ (0. 5*SIN( I 3. 14159/ WAVE ) * IXG-t XSTRT + I 1 .5*WAVE ))) ))> 
ESCWS=ESCWS+ENS 
IF(RSIXNIK) I.LE.G.O ). GC TO 19 
ESCU-ESCD+ENS 
19 CONTiNUE 

I F ( i. MNIN.EQ.O) WRITE(6,3149) 

1 F ( INININ. EQ.O) GO TO 23 ' 

ABESCD=CA13S (ESCb) 

IFIABESCD.LT. 1.0 E-05 ) GO TO 58 

ANES CD=5 7.29 5 7b*ATAN2 I AIMAGI ESCD) * REAL ( ESCD ) ) 

58 CONTINUE 

IFIABESCD .LT. 1.0E-05) ANESCD=0.0 
ESMAGS=CABS< ESCWS ) 

E SANGS= AT AN2 I AIMAGI ESCWS) , REAL IESCWS) ) * 1 80 . 0/ 3 . 14 1 5 92 7 
3149 FORMAT! • %'NO SCATTERED E FIELD WITH SHADOWING*) 

I F C IMNIN.NE.O ) WRITE16. 776 ) E SNAGS , E SAN GS 
776 FORMAT!' * ,'MAG OF SCATT. E FIELD WITH S H AOOW IN G= • , E 15 . 8 , ' P HA SOR 
2AN GLE=* ,615.8) 

IF! IMNIN.NE.O ) WRITE16,2118 ) AB E SC D , AN ESCD 
2118 FORMAT!* *,'SCAT FIELD WITH SHAD. CCNCAVE DOWN ONL Y = * , E 1 5 . 8 , " 

2* PHASOR ANGLE=* ,615.8) 

EWSt I J ) =ESMAGS 
EWCSI I J) =ABESCD 
AW S I IJ ) = ESANGS 
AwCSI IJ)=ANESCD 
GO TO 23 

372 WRITE ( fc » 3 1 52 ) T H ID , THSD 

3152 FORMAT!' NO SPECULAR POINTS FOR TH 1 0= • , El 5 • 8, • AND THSD= * , E 1 5. 8) 
23 WRITE(6,779) 

WRITE <6,779) 

779 FORMAT ( 1 H ) 

17 CONTINUE 


536 


537 


538 


FOR THE PLOTS 
DO 536 I K0= 1 t N AN I 
I ND= I KO-1 
THSO= AOSI IKO ) 

Y ( 1 )=ESNS{ IKO ) 

CALL PLOTITHSD ,Y ,1 ,IND, 50 .0,0.0) 
DO 537 IKQ=1,NANI 
IND= I KO- 1 
TH SD = AOS I IKO) 

Y { 1 ) = E WS ( IKO ) 

CALL PLCT ITHSD ,Y, 1 ,IND, 50 .0,0.0 ) 
DO 538 I K0=1 ,NAN I 
IND= I KO-1 
THSD=AOS( IKO ) 

VI 1 ) = ECDNS( IKO) 

CALL PL OT ITHS D , Y , 1 , I ND , 50 • 0 , 0 .0 ) 
DO 539 I K0=1 , NANI 
I ND= I KO- 1 
THSD= AOS I IKO ) 

YU) -EWCSI IKO) 
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539 


CALL PL0TITHSD,Y,1,IND, 50. 0,0.0) 

00 936 KKRL=1,NANI , 

ANG0i,=FL0ATIKKRL)/2.0 
IPIESNSIKKRLJ.LE. 0.0001) GO TO 936 
DBNS-20 .C*AL0G10( ESNS(KKRD) 

WRITEl 6,937) OBNS, ANGOS 
937 FORMAT l* UBNS= ' ,6 15 . 8, 1 ANGOS=* ,E15.8) 
936 CONTINUE 

00736 KKRL=1,NANI 
ANGQS=FL0AT(KKRL)/2.0 
I F ( IEWSIKKRL) .LE.C.0001) GO TO 736 
DBS=20.0*ALUG10I EWSI KKRL)) 

WRIT EI-6 ,737 ) CBS,ANGOS 

737 FORMAT!* 08S=* , E15.8 , ' ANGOS= » , E15 .8 ) 

736 CONTINUE 
STOP 
END 


FUNCTION RS ( X ) 

COMMON C A, CB , CKA , CKB , PHA , PHB , CC ,C KC , PHC 
C THIS GIVES THE RADIUS OF CURVATURE AT X 

HP=(CA*CKA*COS( <CKA*X) + PHA) ) + (CB*CKB*CO$ I I CKB*X)+PHB) J 
2+lCC*CKC*C0SI ICKC*X)+PHC) ) 

HPP=-( ICA*CKA*CKA*S INI ICKA*X)+PHA) ) + (CB*CKB*CKB*S INI I CKB#X ) +PHB ) ) 
2+ICC*CKC*CKC*SIN( I CKC*X) + PHC ) ) ) . 

RS = UI.O + t HP*HPJ |**1.5)/I-HPPi 
RETURN 
END 

FUNCTION TH(X) 

COMMON CA,CB,CKA»CKB, PHA, PHB, CC,CKC,PhC 

TH=ATAN2I IC A*CKA*COS I (CKA*X) +PHA) ) + ( C B*CKB*COS I I CK8*X J+PHB) ) 
2+(CC*CKC*C0S( (CKC*X) + PHC) ) ,1.0) 

THIS FUNCTION GIVES THE ANGLE BET. THE TANGENT TO HIX) AND THE 
HORIZONTAL 
RETURN 
END 


FUNCTION HIX) 

COHMOR CA,C&,CKA*CKB ,PHA , PH8 , CC ,CKC , PHC 

H= I CA*S INK CKA*X J+PHA I l+(CB*S IN 1 1 CK8*X I +PHB ) ) +CC*SI N 1 1 CKC*X) +PHC ) 

RETURN 

END 



DIMENSION YUO) ,ESSS1360) 

C THIS IS THE TM CASE . 

C THIS PROGRAM USES PHYSICAL OPTICS TO CALCULATE THE BACKSCATTER ING 

C FROM A SEA SURFACE BY DIVIDING SURFACE INTO LIT AND UNLIT REGIONS 

C IN THE LIT REGIONS THE SURFACE CURRENT IN 2NXH 

C GAUSSIAN INTEGRATION USED 

C FOR THIS PROGRAM TO GIVE USEFUL RESULTS THE SURFACE MUST HAVE 

C RADII OF CURVATURE NO LESS THAN 1*WE 

C NSP IS THE HUMBER OF SHADOW POINTS' 

C SURFACE IS DESCRIBED BY AONE*S I N<CCNE*X+PQNE ) +ATWO*SIN( CTWO*X 

C -frPTWO )+ATRE*SIN( CTRE*X+P TRE ) 

C SURFACE UNDER CONSIDERATION LIES BETWEEN ALEP AND REP 
C SN IS THE STEP SIZE TAKEN TO DETERMINE SHADOWING 

C IT MUST BE SMALLER THAN ANY SURFACE FEATURES AND MUST ALSO 

C ALLCW THE LOCATIONOF THE END POINTS OF INTEGRATION WITHIN 

C A SMALL FRACTION OF A WAVELENGTH 

C NANI IS THE NUMBER OF ANGLES (SCATTERING) TO BE EXAMINED 

C MAKE DIMENSIONS OF ESSS , SCANG,EFPA SMALL AS POSSIBLE TO AVDID 

C LAGE U OF CARDS RETURNED 

C NANI SHOULD BE THE DIMENSION OF £SSS , SCANG»EFPA 

NAMELIST/R0N/A8, ANGtDTHS 
DIMENSION SCANG1360) ,EFPA(360) 

COMPLEX S,BINT 

C SCATTER SHADOWING HAS NOT BEEN ACCOUNTED FOR 

COMMON /DOG/AONEt CONE, PONE, AT WO, CT WO ,PTWO, ATRE t CTRE , PTRE 
COMMON /HOG/ G»FHI,THSiWE 
COMMON/ PIG/ SECTORtDX.REPtSECDIO 

COMMON/ GSNN/ GW1 » GW2 , GW3 » GW4, GW5 » GU 1, GU2 , GU3, GU4, GU5 
WE=25.C 

C WE IS THE ELECTRICAL WAVELENGTH 

G=2.0*3. 1415927/WE 
SRTW£=SQRT(WE) 

CX=WE/1 5.0 
A0NE=50.0 

C0NE=2.0*3. 141 5927/800.0 

PCNE=3. 14159/2.0 

ATh‘0-0 « 0 

CTW0=0.0 

PTWD=0.0 

ATRE=0.0 

CTRE=0.0 

PTRE=C.O 

NAN 1 = 360 

SECTCR=WE/2.C 

SECD 10= SEC TOR/ 10*0 

C CONSTANTS FOR GAUSSIAN INTEGRATION 

GW1=0. 2369268 
GW2=0. 47862867 
GW3=0. 568889 
GH4=GW2 
GW5=GW1 

GU1=~C. 9061798 
GU2=-0. 53846931 
GU3=G . 0 
GU4=-GU2 
GU5=-GU1 

C THE ANGLE OF INCIDENCE SHOULD NOT BE GREATER THAN 90 DEG 
TH 1=60. 0*3. 1415927/1 80.0 

C IF THE INCIDENCE ANGLE IS WITHIN TEN DEGREES OF 90 NO SHADOV/ING 

C TAKEN INTO ACCOUNT 

IF(ABS(TFI-1 .5707) .LT. 0.175) GO TO 563 
TANTHI=TAN(THI ) 

DTHI=180.0*T HI/3.141 5 927 
WR l TE ( 6 » 1071 ) DTHI 

1071 FORMAT ( 1 *,* ANG OF INC FROM POS X AXIS =*,E15.8) 



REP=200.0 
ALfcP=-REP 
SN=WE/10.0 
NSP= 1 

DIMENSION SX(IOOO) 

IFIDH(REP) .GT.TANTHI ) GO TO 106 
SX(NSP) =REP 
GO TO 105 
106 SLOPE^TANTHI 

B s: H( REP )-( SLCPE*REP> 

X"REP 

109 X- X- SN 

IF((SLOPE*X)+B.GT.H(XJ 1 GO TO 109 
IP(X .LE.ALEP J GO TO 1000 
$X1NSP)-=X~( SN/2.0J 
105 CONTINUE 

C THIS ABOVE TAKES CARE OF THE FIRST RIGHTENDPOINT 

15 X= SX { NSP ) 

22 X=X-SN 
XN=X-SN 

IFUDMXl.LT.TANTHI ).AND.(DH(XN1 .GT.TANTHIll GO TO 53 
IFCX.GT .ALEP) GO TO 22 
GO TO 92 
53 NSP=NSP+ 1 
SX ( NSP ) =XN 
SLOPE = T ANTH I 

B = H ( S X I N’SP) )-( SLOPE *SX ( NSP ) ) 

X=SX(NSPJ-SN 
29 X=X-SN 

I F € (SLOPE*X)+B.LT.H(X)) GO TO 39 
IF(X.GT.ALEP) GO TO 29 

GO TO 92 
39 NSP=NSP+1 

SX(NSP)=X-(SN/2.0) 

GO TQ 15 
92 NSP=NSP+1 

SX ( NSP I =ALEP 
GO TO 564 
563 sxm=fiep 
SX ( 2 ) =AL EP 
NSP = 2 

564 CONTINUE 

C LAST VALUE IN SX(J) IS ALEP 

WRITE (6*101 » ( K,SX(K) , K=1 f NSP) 

101 FORMAT ( * 1 SXl* , 14, » ) = * ,E15.8) 

DO 317 JNX = 1 , NAN I 

THS = FLOA T(JNX>*(0. 8726646 E-021 

DTHS=180 .O+THS/3. 1415927 

SCANGI JNX)=DTHS 

S=CMPLX( 0.0,0.01 

KKN= 1 

10 CONTINUE 

ALCW=SX( KKN+1) 

AUPP=SXI KKN ) 

S=S+BINT( ALOW,AUPP) 

KKN=KKN+2 ' 

IF { (KKN.LT.NSPl .AND. ( (KKN+1 1 .LT.NSPI I GOTO 10 
C TO CONVERT TO TRUE SCATTERED E FIELD FOR EINC OF UNITY HAG, 

S=CMPLX ( -0.7071 1 ,-0.707111 *S/SRT WE 
AB=CABS( S) 

ESSS ( JNX )=A8 

ANG=180.0*ATAN2( A l MAG (S) »REAL ( S ) 1/3.1415927 
EFPAt JNX )=ANG 
317 CONTINUE 

DO 531 J K= 1 » NAN I 


93 
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^E=ESSS( JK) 

DB=20.0*ALOG10(E ) 

A=EFPM JK) 

AS=SCANG( JK) 

531 WRITE < 6 » 532 I AS,E , A »DB 

532 FORMAT ( ' SCAT ANG FROM H0RIZ=' ,E15.8»* MAG OF E FIELD=», 

2 E15.B, 1 PHASE ANG= * , E15 .8 , * 0B=*,E15.8) 

DO 535 I KE=1 »NAN I 
INO=IKE-i 

THSD=FLOAT( IKE 1/2.00 
Y ( 1 ) =ESSS( IKE ) 

535 CALL PLOT < THSD» Y, 1 , IND, 50 .0,0 .0 ) 

GO TO 1002 
1000 WR I TE ( 6, 1592 ) 

1592 FORMAT (• SURFACE IS NOT . I LLUM 1NAT ED • ) 

1002 CGNTINUE 
STOP 
END 


COMMON / DOG/ AONE , CONE » PONE » AT WO » CT WO » PT WO , AT RE »C TRE , PTRE 

H-AONE*S INI CONE* X + PONE ) *ATWO*S IN (CTWO*X+PTWO l+ATRE*S INICTRE* X+PTRE 

2) 

RETURN 

END 


FUNCTION DH IX) 

COMMON /DCG/AGNE ,CONE,PONE, AT WO, CT WO , PT WO, ATRE ,C TRE, PTRE 
DH=AONE*CONE*COS ( CONE*X+ PONE ) +AT WO *CTWO*COS ( C TWO *X+P TWO ) 
2 + ATRE*CTRE*COS(CTRE*X-»-PTRE) 

RETURN 

END 


FUNCTION B I NT ( XX , Y Y ) 

XX IS LOWER LIMIT OF INTEGRATION, YY IS UPPER LIMIT 
PHYSICAL OPTICS RADIATION INTEGRAL WITH PLANE WAVE INCIDENT 
TM CASE 

COMPLEX S » BI NT 
COMPLEX GASS5 

COMMON /HOG/ GtTHIfTHSfWE 
COMMON/PIG/ SECTOR ,DX,REP,SECD 10 

BRE^fC INTEGRAL FROM XX TO YY INTO SHALLER SEGMENTS - OF LENGTH 
SECTOR AND INTEGRATE OVER EACH SEGMENT USING GAUSSIAN INTEGRATION 
S=CKPLXIO. 0,0.0) 

LDS=INT( I YY-XX I / SECTOR ) 

IF (LOS. EQ.O ) GO TO 10 
DO *00 I NJ= 1 , LDS 
UL=XX+( FLOAT ( INJ)*SECTOR) 

ALL=*-XX *( FLOAT (IN J-l )*SEC TORT 
100 S=S+GAS$5( ALL,UL) 

C NOW TO GET LAST FRACTION OF SEGMENT LEFT OVER FROM SURFACE SEGMENTATION 

S=S+GAS S5IXX+C FLOAT! LDS ) *S ECTOR) ,YY) 

GO TO 50 

10 S=GASS5 I XX , YY ) 

50 CONTINUE 
BINT=S 
RETURN 
END 
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FUNCTION GASS 5 (XL, XU) 

COMPLEX GASS5 »FTBI 

FIFTH ORDER GAUSSIN INTEGRATION 

XL IS LOWER LIMIT, XU IS UPPER LIMIT 

XU-XL IS LESS THAN OR EQUAL TO SECTOR 

COMMON/ GSNN/GWl , GW2 , GW3 ,GW4 , GW5 , GU X , GU2 , GU3 , GU4, GU5 

DV0FEP=(XU-XL)/2.0 

0VSMEP=(XU+XL)/2.0 

XU5=GU5*DVOF£P+DVSMEP 

XU4=GU4«0VDFEP+0VSMEP 

XU3=GU3*DVDFEP+DVSMEP 

XU2=GU2*0V0FEP+DVSMEP 

XU1=GU1*OVOFEP+DVSMEP 

GASS5=DV0FEP*(GWX*FTBI( XUl )+GW2*FTBI { XU2 ) +GW3*FTBI (XU3) 
2 +GW4*FTBI(XU4H-GH5*FT8I (XU5) ) 

RETURN 

END 


FUNCTION FTBIIX) 

COMPLEX FTBl 

THIS IS THE FUNCTLGN TO BE INTEGRATED 
THIS IS FOR THE TM CASE 
COMMON/HGG/G , THI ,THS,WE 
COfMON/PIG/ SECTOR ,DX , REP, SECDIO * 

GCC^G^-I COS (THI ) +CQS ( THS ) ) 

GSS=G*(SIN(THI)+SIN<THS> ) 

RCK=REP-(2.0*WE) 

FTBI=SIN(THI-ATAN(DH(X) ) )*SQRTll.O + ( CH(X>**2 ))* 

2 CEXP<CMPLX(0.0,( (X*GCC>+(H<X)*GSS) )) ) 

C THE FOLLOWING ACCOUNTS FOR TAPERING 
ABSX=ABS ( X ) 

IFIA8SX-RCK) 1500,1500,2000 
2000 IFtX.LE. (WE-REP) ) FT8 l=CMPLX ( 0 .0 , 0 .0 ) 

IFCX.GE. (REP-WE) ) FTBI=CMPLX (0.0,0 .0 ) 

IF((X.GT.( WE-REP) ) . AND. ( X, LE. ( ( 2 .0*WE)-REP ) > ) 

2 FTBI=FTBI*(0.5+( 0.5*SIN( (G/2.0) $( X- ( ( 1 . 5*WE )-REP) ) ) )) 
IF( ( X.LT. (REP-WE )) .AND. (X.GT.( REP- (2. 0*WE) )) ) 

2 FTB!=FTBI$(0.5-(0.5*SIN( ( G/2.0) *( X- (REP— ( 1.5*WE )) ) ) ) ) 
1500 CONTINUE 
RETURN 
END 
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C THIS IS THE TE CASE 

C THIS PROGRAM USES PHYSICAL OPTICS TO CALCULATE THE SACKSCATTERING 

C FROM A SEA SURFACE BY DIVIDING SURFACE INTO LIT AND UNLIT REGIONS 

C IN THE LIT REGIONS THE SURFACE CURRENT IN 2NXH 

C GAUSSIAN INTEGRATION USED 

C FOR THIS PROGRAM TO GIVE USEFUL RESULTS THE SURFACE MUST HAVE 

C RADII OF CURVATURE NO LESS THAN L*WE 

C NSP IS THE NUMBER OF SHADOW POINTS 

C SURFACE IS DESCRIBED BY AONE*SIN( CONE*X+PONE ) -«-ATWO* SIN! CTWO*X 

C +PTWO)+ATRE*SINICTRE*X*PTRE) 

C SURFACE UNDER CONSIDERATION LIES BETWEEN ALEP AND REP 

C SN IS THE STEP SIZE TAKEN TO DETERMINE SHADOWING 
C IT MUST BE SMALLER THAN ANY SURFACE FEATURES AND MUST ALSO 

C ALLOW THE LOCATIONOF THE END POINTS OF INTEGRATION WITHIN 

C A SMALL FRACTION OF A WAVELENGTH 

.C NANI IS THE NUMBER OF ANGLES SCATTERING) TO BE EXAMINED 

C MAKE DIMENSIONS OF ESSS , SCANG»EFPA SMALL AS POSSIBLE TO AVOID 
C LAGE # OF CARDS RETURNED 

C NANI SHOULD BE THE DIMENSION OF ESSS ,SCANG,EFPA 

DIMENSION YUO), ESSSI360) 

NAMEt I ST/RON/AB ,ANG ,DTHS 
DIMENSION SCANGC360) ,EFPAI 360) 

COMPLEX S t B I NT . 

C SCATTER SHADOWING HAS NOT BEEN ACCOUNTED FOR 

COMMON /DOG/AONE, CONE, PONE, ATWO,CTWO, PTWO, ATRE,CTRE,PTRE 
COMMON /HOG/ G,THI,THS,WE 
COMMON/PIG/ SECTOR, OX, REP, SECD10 

COMMON/ GSNN/ GW I, GW2t GW 3, GW4, GW5»GUL»GU2,GU3,GU4,GU5 
WE=25.0 

C WE IS THE ELECTRICAL WAVELENGTH 

G=2 .0*3. 1415927/WE 
SRTWE=SQRTl HE) 

DX=WE/15.0 
AON E— 40 .0 

C0NE = 2. 0*3. 1415927/200.0 

PONE-O.O 

ATW0=0.0 

CTW0=0. 0 

PTWO=0 .0 

ATRE=0 .0 

CTRF=0. 0 

PTRE=0 .0 

NANI =360 

$ECTQR=WE/2.0 

SECD10=SECTOR/10 .0 

C CONSTANTS FOR GAUSSIAN INTEGRATION 

GWI=0. 2369268 
GW2 =0.47862 867 
GW3 =0.568889 
GW4=GW2 
GW5=GW1 

GUI =-0.9061 798 
GU2=-0. 53846931 
GU3=0 .0 
GU4=-GU2 
GU5=-GU1 

C THE ANGLE OF INCIDENCE SHOULD NOT BE GREATER THAN 90 DEG 
THI=60. 0*3. 141 5927/ 180.C 

IF THE INCIDENCE ANGLE IS WITHIN TEN DEGREES OF 90 NO SHADDOWING 
TAKEN INTO ACCOUNT 

IF(ABS(THI-1. 5707). LT. 0.175) GO TO 563 
TANTHI=TAN( THI 1 
DTHI=180.0*TH 1/3. 1415927 
WRITE (6, 1071) DTHI 
1071 FORMAT ( * », * ANG OF INC FROM POS XAXIS =*,£15.81 
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REP=200. 0 
AL EP=-REP 
SN=W E/ 10 .0 
NSP=1 

DIMENSION SXI 1000) 

I F I DH I REP } . GT ,T ANTH I I GO TO 106 
SXINSP)=REP 
GO TO 105 
i 06 SLOPE =T ANTH I 

B=HI REP)-< SLOPE+REP) 

X = REP 

109 X=X-SN 

IFK SLOPE*X)+B.GT.H( X) I GO TO 109 
1F(X .LE.ALEPI GO TO 1000 
SX( NSP) =X-ISN/2.0 ) 

105 CONTINUE 

C THIS ABOVE TAKES CARE OF THE FIRST R IGHTENDPOI NT 
15 X^SXINSP) 

22 X=X-SN 
XN=X-SN 

IF KDH(X) .LT.TANTHM . AND, (OH(XN) .GT .TANTHI ) I GO TO 53 
IF f X.GT.ALEP) GO TO 22 
GO TO 92 
53 NSP=NSP+1 
SXINSP ) =XN 
S LOPE=T ANT H I 

B=HISX(NSPI )-( SLOPE* SX (NSP I ) 

X* SXI NSP ) *-SN 
29 X=X-SN 

IF I ( SLOPE+X I + B. LT . H IX ) ) GO TO 39 
I F ( X.GT.ALEP! GO TO 29 
GO TO 92 
? NS P= NS P+1 

SXINSP) = X-( SN/2 .0) 

GO TO 15 
92 NSP=NSP+1 

SXtNSP)=ALEP 
GO TO 564 
563 SX II ) =REP 
SXI 2) =ALE P 
NSP = 2 

564 CONTINUE 

SURFACE IS NOW SEPERATEO INTO LIT AND UNLIT ZONES 
LAST VALUE IN SX(J) IS ALEP 
WRITE 16,101) (KtSX (KI,K=l ,NSP) 

101 FORMAT! • » , • SXI • , 1 4 , • ) *• ,E 1 5. 8) 

C THE FOLLOWING F IN DS THE SCATTERED FIELDS DUE TO THE LIT ZONES 

DO 31 7 JNX=1 ,NANI 
T HS 3 FLO ATI JN X ) *1 0 • 87 26646 E-02) 

DTHS=180.0*THS/ 3 .1415927 
SCANGI JNX) =DTHS 
S*CMPLX(0. 0,0.0) 

KKN=1 

10 CONTINUE 

ALOW-SX I KKN + 1 1 
AUPP=SX!KKN) 

S=S+B I NT I ALOW , AUPP) 

KKN=KKN+2 

IF ( t KKN. LT • NSP) • AND. I IKKN+l ) .LT .NSP ) ) GO TO 10 
C TO CONVERT TO TRUE SCATTERED H FIELD FOR HINC OF UNITY MAG 

S=S*CMPLX(0. 7071:1, 0.70711 l/SRTWE 
AB=CA B S I S) 

DB=2C .0*ALOGIOIAB) 

ESSSI JNX l*AB 

ANG=180. 0*ATAN2I AIMAGIS) ,RE AL IS ) 1/3.1415927 
EFPAI JNX )=ANG 

WRITE (6, 143) OTHS,AB» ANG, DB 
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143 FORMA T( * SCATTERING ANG*,E15.8** HAG=* , E15 .8 * » PHASE AN GLE=*, 
2 El 5 .8 , 1 DB=',E15.8> 

317 CONTINUE 

DO 531 JK=l,NANI 
E=ESSSIJK» 

A=EFPA( JK) 

AS-SCANG( JK I 

531 WRITE (6,5321 A$,E,A 

532 FORMAT ( ' • ,» SCAT ANG FROM HORIZ=* ,£15.8,* MAG OF H FIELD**, 

2 E15.8,' PHASE ANG=*,E15.8) 

DO 535 IKE=1,NANI 
I ND = I KE- 1 

THSD=FLOAT( IKE) /2.00 
Y( l l=ESSSUKE) 

535 CALL PLOT < THS 0 » Y , 1 , I ND, 50 .0 , 0 .0 ) 

GO TO 1002 
1000 WRITE (6,1592 > 

1592 FORMAT( • SURFACE IS NOT I LLUMI NAT ED* 1 
1002 CONTINUE 
STOP 
END 


rnMMDN /DOG /AONE , CONE , PONE, ATWO,CTVJO , PTWO ,ATRE ,CTRE , PT RE r. TOC 

H= AONE*S^N^ CON £*X +PON E ) +ATW0 *S IN ( CTWQ*X*P TWO >+ATRE*SIN(C TRE* X+PTRE 

2 > 

RETURN 

END 


Z 4-ATR£*CTRE*COS (CTRE*X+PTRE) 

RETURN 

END 


FUNCTION B INT( XX, YY) 

XX IS LOWER LIMIT OF INT EGR AT ION, YY IS UPPER LIMIT 
PHYSICAL OPTICS RADIATION INTEGRAL WITH PLANE WAVE INCIDENT 
TM CASE 

COMPLEX S, BINT 
COMPLEX GASS5 
COMMON /HOG/ G» THI ,THS» WE 
COMMON/PIG/ SECTOR, OX, REP, SECD10 

BREAK INTEGRAL FROM XX TO YY INTO SMALLER SEGMENTS OF LENGTH 
SECTOR AND INTEGRATE OVER EACH SEGMENT USING GAUSSIAN INTEGRATION 
S=CMPLXt 0.0,0.01 
LDS=INT( ( YY- XX) /SEC TORI 
IFILDS. EQ.O I GO TO 10 
DO 100 INJ*1,L0S 
UL*XX + (FLOAT< INJ ) +SECTOR) 

ALL=XX+ < FLOAT (I NJ-1 )*S ECTOR I 
100 S=S*GA$S5(ALL,UL1 

S=S*GASS5(XX+( FLOAT (LDS>*SECTOR ) , YY» 

GO TO 50 

10 S=GA$S5< XX, YY) 

50 CONTINUE 
BINT=S 
RETURN 
END 
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FUNCTION GASS5 (XL, XU I 

COMPLEX GASS5 ,FT8 I 

FIFTH ORDER GAUSSIN INTEGRATION 

XL IS LOWER LIMIT »XU IS UPPER LIMIT 

XU- XL IS LESS THAN OR EQUAL TO SECTOR 

COMMON/ GSNN/GW1, GWZ, GW 3, GW«,GW5,GU1.,GU2,GU3,GU4,GU5 

OVDFEP=( XU-XLI /2.0 

DVSMEP= ( XU* XL) /2. 0 

XU5=GU5*DVDFEP*DVStf EP 

XU4=GU4*DV0FEP+DVSMEP 

XU3=GU3*DVDFEP*DVSMEP 

XU2=GU2*DVDFEP+OVSMEP 

XUl=GUI*DVOFEP*OVSMEP 

GASS5=DVOFE P* (GW I*F T BI ( XU i) * GW2*F T8 1 ( XU2I+G W3*F TBI < XU3I 
2 *GW<»*FTBI (XU4) *GW5+FTBI (XU5JI 
RETURN 
END 


FUNCTION FT8IIX1 
COMPLEX FTBI 

THIS IS THE FUNCTION TO BE INTEGRATED 
THIS IS FOR THE TM CASE 
COMMON/HOG/G, THI , THS , WE 
COMMON/PIG/ S ECTOR, DX, REP, SECDIO 
GCC =G*( COSf THI )+COS( THS) ) 

GSS = G*(SIN(THI >+SIN(THS) ) 

RCK=REP-(2.0*WE) 

FTBI*SIN( THS-ATAN(DH( X) ) )*SQRTI I . 0+ ( DH ( X ) **2 ) )* 

2 CEXP(CMPLX(0,0, I ( X *GCC ) +( HI X)*G$S) 1 )) 

C THE FOLLOWING ACCOUNTS FOR TAPP£RING 

ABS X= AB S( X ) 

l F ( ABSX -RCK ) 1500,1500,2000 
2000 IF( X. LE. ( WE-REP) I FT0 I =CMPLX (0.0, 0.0 J 
IF(X .GE.(REP-WEI) F TB l =CMPL X( 0. 0, 0. 0> 

IF ( (X.GT. (WE-REP) ) . AND. ( X . L E . ( (< .0*WE)-REP) ) I 
2 FTB I=FT8I* ( 0. 5* ( 0. 5*S I N( (G/2 .0* * ( X- ( (I . 5*WE)-REP| III) 
IF{ (X .LT .IREP-WE) ) . AN D . ( X .GT . ( REP- ( 2 . 0*WE > ) ) ) 

2 FTBI =FTBI*(0.5-(0.5+S INI (G/2 .0 ) * ( X- { REP- ( 1 ,5*W E 1) ))) ) . 
1500 CONTINUE 
RETURN 
END 
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THIS IS A METHOD OF MOMENTS SOLUTION 
TM POLARIZATION SYMMETRIC MATRIX 
N SUB SEGMENTS HAV0 N MIDPOINTS 

NSUB.IS THE SUBSCRIPT WHICH COUNTS THE END POINTS 
N IS THE SUBSCRIPT WHICH COUNTS THE MIDPOINTS 
WATCH MAX SLOPE SO THAT THE X INCREMENTS ARE SMALL ENOUGH 
THE REGION UNDER CONSIDERATION LIES BETWEEN -EP AND EP 
DIMENSION Y{ 10) ,CMC(360) 

CGMPLEX SNN.SST 
COMPLEX FSS 

DOUBLE PRECISION DAL , DDX , DDC2 , DDC , DALC , DR 
COMPLEX FINCI30 ),5TS 

COMMON /PIG/ AONE,CONE f PCNE,ATWO,CTWG, PTWOiN 
COMPLEX A HAN 20 

COMPLEX F(300),SIA515O)»SS,T 
COMPLEX FIN 
DIMENSION X ( 300 ) 

DIMENSION XM I D (300 ) 

COMPLEX STO 

WE IS THE ELECTRICAL WAVELENGTH 
WE=25.0 

G=6-2831853 /WE 

S I S = SQRT (WE)*CMPLX(l .0 , 1 . 0 » * ( +0 . 70 7 10 7 ) / 3. 14 1 592 7 
DC=NE/10.0 
DX = DC / 1 OC'Oi • 0 
DC2 -DC/2.0 
EP=200. 0 
AP I =3 • 1 A 1 592 7 

THE FOLLOWING CONSTANTS DEFINE THE SURFACE 
A0NE-25.0 

CCNE=2 .0*3. 1415927/20C.0 
PON E=G • 0 
AT^Q=G .0 
CTWO=C .0 
P T WO = 0 • C 
CALL SCL0K1 

THE FOLLOWING BREAKS THE SURFACE INTO SEGMENTS DC CENTIMETERS LONG 
BY LINE INTEGRATION USING STEPS OF LENGTH DX FOR THF INTEGRATION 
NSU8=I 
X ( NS UB ) =*-E P 
DOC = DBLE< DC) 

DDX= CBL E (DX) 

DDC2=0BL E( DC2 ) 

1002 DAL=0.000 00 

dr=cele(x(nsub) ) 

1001 DR=DR+DDX 
R=SNGL ( DR ) 

DALQ=DAL 

C AL = CAL+ { CDX + DSQRT t 1 .00 00 + ( ( DB L E ( QH( R ) ) ) ** 2 ) ) ) 

IF( ((DDC2-DAL) .Lt.C.OO CO ) . AND. ( ( DDC2-CAL0 ) . GE.O .00 00)) 

2 X M I D ( NSUB ) =K 
IF( DAL.LT .CDCJGQ TO 1001 
NSUB=NSU9+1 
X(NSU8J=R 
AL=SNGL( CAL ) 

WRITE (6,352) A L » NSUB 
352 FORMAT ( • * , ■ AL = • , E 1 5. 8 , » NSUB=»,IA) 

l F (R.LT .EP ) GO TO 1002 
T I ME = RC LCK I ( l.C) 

WR IT E ( 6 , 32 76 ) TIME 

3276 FORMAT ( • * , • T I ME = * » F 10 . 6 , • S ECCNDS » ) 

N=N$UB-1 

DO 100 <=f J=1,NSUB 
IF (J.EQ.NSUB) XMID(NSUB) =0 .0 
XXX=X( J ) 

XMO~XMIO( J ) 


TOO 
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1004 WRITE ( 6 t 1003 > XXX,XMD,J 

1003 FORMAT ( 6h X ( J ) = , E 1 5 . 8 , 9H XMI D ( J ) = , E 1 5 . 8 , 3H J=,I3) 

C THIS ENDS THE SURFACE SUBDIVISION 

NMO-N- 1 
NM3=N-3 ’ 

C PI MENS I ON OF S IS NIN+D/2 

C OIMeHSlON OF F INC* F IS N 

0 P I F =0 .7853982* 

EE = 2, 71828 
GA=G*DC/ ( 2.0*EE > 

C SNN IS THE DIAGONAL ELEMENT OF THE INPUT MATRIX 
SNN=AHAN20(GA) 

WRJTE (6,400) SNN 
400 FORMAT ( 5H SNN=,2E15.8J 
DO 100 N J= 1 »N 
NJPD=NJ+I 

S( ISUB(N4,NJ) )=SNN 

C THIS FINDS ELEMENTS ON THE DIAGONAL 

IF (NJPO.GT.tf) GO TO 100 
DO 100 na=njpo,n 

C THIS FINDS OFF DIAGONAL ELEMENTS 

XM=XMID ( N J ) 

XN=XM1 D < NA) 

RHO=SQRT ( ( ( XN-XM )**2) + ( ( H < XN ) - H( XM ) ) ** 2 ) ) 

RHG=RHO*G 

SC I SLB ( N J , NA ) ) = AHAN20 ( RHG ) 

ICO CONTINUE 

C THIS COMPLETES THE FILLlN OF THE MATRIX 

C THIS BEGINS THE CONVERSION TO UPPER TRIANGULAR MATRIX 

S( 1) = CSQRT( S< 1 ) ) 

DO 1 K = 2 , N 
1 S(K)=S(KI/S(1) 

DO l 1 = 2 , N 
IMG= I —1 
IPO=I+I 

T=CMPLX( 0.0, 0.0) 

DO 3 L = 1 , I MO 

LI=<L*N>-< ( ( <L-U*L)/2)+N-I) 

3 T = T+{S( LI ) **2 ) 

ii = n*N»-((((i-i)*n/2 j+n-i > 

St 11) =C SQR T I S( II )-T) 

1 F ( 1P0.GT.N ) GOTO 2 
DO 5 J = I PO , N 
T=CMPLXI0.0 ,0*0) 

DC 6 M= 1 , IMU 

MI = (M*N)-( ( ( (M-l )*M)/2 )+N~I > 

MJ=(M*N)- { ( (MMM-l) ) /2) +N-J) 

6 T=T*(S(MJI*S(MI ) ) 

IJ = ( I*N)-( ( ( ( I -1 J * I ) / 2 ) + N- J ) 

S( I J ) = ( S( 1 J )- T )/S( I I ) 

CONTINUE 

C THIS ENDS THE CCNVERSICN TO UPPER TRIANGULAR MATRIX 

WRITE (6,1222) N,WE 
1222 FOWTOh M=»I3»4H WE=,E15.8) 

TH=60. 0*3. 1415927/180. C 
THXXD= 1 80. O+TH/3. 1415927 
WRIT£ (6,9333) THXXO 
9333 FORMAT ( 9H INC ANG=,E15.8) 

C TH IS THE ANGLE OF INCIDENCE FROM THE HORIZONTAL 

STH=SrIN(TH) 

CTH=COS(TH) 

C THIS FINOS THE INCIDENT FIELD ION THE N J TH SEGMENT 

DO 455 HJ=1 , N 
ENJ=FLOAT(NJ ) 

XM=XMID(NJ> 

F(NJ)=CEXP(CMPLX(0.0,G*( ( XM*CTH ) + ( H ( X 1 ) *ST H) ) ) ) 
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TAPERED ILLUMINATION ***************************** *********** 


IFCXM.LE.t (WEM.C >-EP) > F ( N J > =CMPLX ( 0 . 0, 0 .0 ) 

I F( XM.GE.I EP-( 1 .0*WE ) » ) F ( N J )=CMPLX (0.0 ,0.0 ) 

IF ( (XM.GT. ( ( l .C*WE )-EP ) ) .AND.I XM .L E . < { 2. 0*WE )-EP ) )) 

2 F( NJ)=F(\J )*{ 0.5MO.5*SIN( (G/2.0 )*(XM -{ ( 1.5*WE >-EP > 1 I ) ) 

IFUxM .GE.(EP-(2.0*WE) M.ANU.IXM . L T. ( EP- ( 1 . 0*WE ) ) )) 

2 F(NJ ) = F( NJ )*(O.5-(0.5*SlNUG/2. D)*( XM - { EP- ( 1 . 5*WE ) ) ) ) ) ) 

455 CONTINUE 

WRITE<6,2948) ( NJ ,F I N J ) , NJ =l',N I 

294$ FURM AT ( * »,» INC FIELD F { • , I 4 , • ) = • f 2 E 1 5. 8 ) 

C THIS BEGINS THE BACK SUBSTUT 1UN 

F< 1>=F< 1 » / S < 1) 

DO 10 1=2, N 

I MOM- l 

T=CNPLX< 0. 0, 0.0) 

DO II L=1 , IMG 

Ll=(L*N)-{ (( ( L- 1 ) * L ) / 2 ) +N“ I ) 

11 TsTHSILI)*FILII 

I I = l I*N ) -C t I I 1-1 1* I ) / 2 ) +N- 1 ) 

10 F(I)=(F< II-TI/SUI ) 

NN= ( N* ( N + 1 ) ) / 2 
F(N)=F(N)/S( NN ) 

NMC=N- I 
DC 25 1=1, NMO 
K=N- 1 
KPO=K+ 1 

TaCMPLX 10 .0,0.0 ) 

DC 26 L=KPG,N 

Kl = (K*N >-(((< K-l )*K) /21+N-L) 

26 T::T + (SIKU*FIL1) 

KVC=< K*N)-( ( ( (K-l )*K) / 21+N-K1 
FtKJ = C FIKJ-T )/S(KK ) 

25 CONTINUE 

C THIS ENDS THE BACK SUBSTITUTIONS 

DO 491 K= 1 , N 
STT = CA8S ( F ( K I 1 
STO=F (K I 

AHNN=ATAN2 (/UMAGtFIK 1 ),REAL (F(K) ) 1*180.0/3.1415927 

491 WRITE(6,^92) K ,STO,STT,ANNN 

492 FORM AT I * • , • F ( » , 14 , • ) = ' , 2E 1 5. 8 , » DR ■ , • AMP = ' , E 1 5 . 8 , • AT ANGLE = ' , 
2 G15.8) 

DO 317 JNX = 1 ,360 

TH=0* 872.6E46E5E-02 *FLOAT(JNX) 

T --CMPLX (0.0,0. 0 1 
□ □ 310 1=1, N 
XN = XM I D ( I 1 

310 T=T+ ( ( F ( I 1 * CEX P ( CMP LX ( C .0, G* ( (XN*COS ( TH 1 ) +{ H( XN )*SIN( TH ) ) 1 ) ))) 

C THIS CORRECTS T TO TRUE SCATTERED FIELD 

T=STS*T 
CM=C ARS ( T ) 

CMC ( JNX 1=CM 

CANG=57 .296* AT AN2 ( A IMAG( T ) , REAL ( T ) ) 

T HP =T hf* 5 7 • 2 96 
DB=20.0*ALOG10(CM) 

317 WRITE (6,312) C M,C ANG , THU, DU 

312 FORMAT (18H RELATIVE E F IEL D= , El 5. 8 , 7H AN GL E= , E 1 5 . 8 , 

2 23H ANGLE FROM HOPI ZUNT AL= , E 1 5 . 8 , 7H DB=,E15.8) 
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DO 576- I KE = 1 » 360 
THSO=FLDAT I I KE)/2.0 
IN0=IKE-1 
Y ( 1 )=C*C ( IKE) 

576 CALL PLOT! THSD,Y, l ,IND*5G .0,0.0) 

STOP 

END 

FUNCTION H(X) 

THIS DEFINES THE SURFACE 

COMMON /PIG/ /'ONE, CONE, PONE, ATWO , CTWO ,PTWO ,N 
H=AGNE*SlN(CUNE*X+PONE I +ATWQ*S IN ( CTWO*X+PTWU) 

RETURN 

END 

FUNCTION DHIX) 

DH(X) IS THE DERIV. OF H(X) 

COMMON /PIG/ ACNE, CONE, PONE, ATWO, CTWO, PTWO, N 
DH=ApNE*CONE*COS ICON E*X+PQN E ) +AT WO *C TWO*CO SI C TWO* X+ PT WO ) 
RETURN 
END 

FUNCTION l SUBt J » K ) 

COMMON /PIG/ AONE, CONE , PONE, AT WO , C TWO , PT WO , N . 

THIS CONVERTS ELEMENTS OF UPPER TRIANGULAR MATRIX 
I SU8 = (N*J) -I II I J-l ) * J ) /2 ) +N-K ) 

ARRAY COUNTING LEFT TO RIGHT STARTING WITH FIRST 
RETURN 
END 


TO A LINEAR 
ROW 
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THIS ISA METHOD OF MOMENTS SOLUTION FOR BI STATIC SCATT TM CASE 
GAUSSIAN INTEGRATION IS USED TO CALCULATE THE MATRIX ELEMENTS 
UNIT INCIDENT ELECTRIC FIELD- IS ASSUMED, OF COURSE THIS IS MODIFIED 
NEAR THE ENDPOINTS OF THE SURFACE BY ILLUMINATION TAPPERING 
NSU8 SEGMENTS HAVE N MIDPOINTS 

NSU3 IS THE SUBSCRIPT WHICH COUNTS THE END POINTS 
N IS THE SUBSCRIPT WHICH COUNTS THE MIDPOINTS 

WATCH MAX SLOPE SO THAT THE X INCREMENTS ARE SMALL ENOUGH 
'THE SURFACE UNDER CONSIDERATION LIES BETWEEN -EP AND +EP 
THE ARRAY XM(J) CONTAINS THE X COORDINATES OF THE MIDPOINTS OF THE 
SEGMENTS, XM( I) IS THE MIDPOINT OF THE I'TH SEGMENT 
THE ARRAY X(J) CONTAINS THE X COORDINATES OF THE ENDPOINTS OF THE 
SURFACE SEGMENTS, X(I),X( 1+1) ARE THE LOWER AND UPPER X COORDINATES 
OF THE ENDPOINTS OF THE I'TH SEGMENT 

PHASE REFFERENCE IS AT THE ORIGIN OF. THE COORDINATE SYSTEM 
COMPLEX SNN,SST 
COMPLEX S 

DIMENSION Y ( 1 0 1 , CMC ( 360 J 

NAMELIST/D/ WE , C P * THXXD , AONE , CONE ,PONE , AT WO , CTWO , PTWO , N 
NAMELIST /E/P,XMIO 
COMPLEX FSS 
COMPLEX STS 

COMMON /PIG / AONE, CONE , PONE, ATWO, CTWO , PTWO, N 
COMPLEX C ( 236 * 236 ) 

COMPLEX F (2361 ,SS,T,CTEST 

THE DIMENSIONS OF C AND F MUST BE COMMENSURATE 

THAT IS C ( L , L ) FILJ 

COMPLEX FIN 
COMPLEX HAN2 
DIMENSION X { 500 } 

DIMENSION X M ( 5 0 0 ) 

C THE FOLLOWING CONSTANTS DESCRIBE THE SURFACE 

ACNE=-50 . 0 
C0NE=6. 28318/800.0 
PGNc = 3 . 1415927/2.0 
A TWO=0 . 0 
C TW0=0 • 0 
P TW0=0 . 0 

C WE IS THE ELECTRICAL WAVELENGTH 

WE=25.C 

G=6. 2831853 /WE 

DC=WE/10.0 

DX=DC/ 1 COO • 0 

0C2=DC/2.0 

EP=20Q . 0 

C THE FOLLOWING BREAKS THE SURFACE INTO SEGMENTS DC CENTIMETERS LONG 

C BY LINE INTEGRATION USING STEPS OF LENGTH DX FQR THE INTEGRATION 

NSU3 = I 
X( NSL3)=-EP 
1002 A L = 0 . 00 0 
R = X(MSUB) 

1001 R=R+DX 
ALO= AL 

AL = AL*( DX*SQRT( 1 ■. 0 + ( DH ( R ) **2 ) J ) 

I F ( < (DC2-AL J ,LE .0 .0 ) . AND. ( ( DC2-ALC) .GT .0. 0) ) XM(NSUB)=R 
IF< AL.LT.DOGO TO 1001 
WRITE ( 6,352) A L , NSUB 

352 FORMAT ( * ' , ' AL= ' , E 15 . 8 , ' NSUB=',I4) 

N SUB=NSUB + 1 
X(NSU3)=R 

IF (K.LT.EP) GO TO 1002 
N=N SUB- 1 

DO 1004 J=1 , NSUB 

IF ( J , EO.NSUB ) XM ( NSUB ) =0 . 0 

XXX=X( J ) 

XMD= XM( J ) 
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100 A WRITE ( 6,1003) XXX,X'MD,J 

1003 FORMflT (6H X ( J ) = , E 1 5 . 8 , 9H XM{ J ) =, E15 . 8, 3H J=,13J 

C THIS ENDS THE SURFACE SUBDIVISION 

NMO=tf-l 
NM3=^-3 

C DIMENSION OF FIN C,F IS N 

DPI F=C • 7853982 
EE = 2. 71828 
GA=G*0C/(2.0*EE ) 

C SHftf IS THE DIAGONAL ELEMENT OF THE INPUT MATRIX 
SNN=HAN2<GA)*0C 
WRITE l 6,400 ) SNN 
400 FORMAT (5H SNN=,2E15.8) 

DO 100 NJ=1,N 
C(NJ,NJ )=SNN 
100 CONTINUE 

C CONSTANTS FOR GAUSSIAN INTEGRATION 5 TH ORDER 

GUl=-0. 9C61798 
GU2=-0. 53846931 
GU3-0.O 
GU4=-GU2 
GU5 =-GU 1 
GW1=0. 2369268 
GW5=0. 2369268 
GW4 =0 .47862867 
GW2=0. 47862667 
GW3=0 • 5688888 
DO 3361 MR= 1 . N 


XMM=XM(rtft) 

HXMM=H( XKMJ 
DO 3361 MC=1,N 
IF (MC.EQ.MR) GO TO 3361 
EPL=X(MC> 

EPU = X ( MC + 1 ) 


DVDFEP= (EPU-EPU /2.0 
DVSMEP=(EPU+EPL)/2.0 
XU5=GU5*DVDFEP+OVSMEP 
XU1=GU1*DVDFEP+DVSMEP 
XU2 = GU2<=DVDFEP + DVSMEP 
XU3=GU3*DVDFEP+DVSMEP 
XU4=GU4*DVUFFP+DVSMEP 
C ( MR , MC ) = DVDFE P* ( 

2+GW1 *HAN2( G*SQRT( ( ( XU 1-XMM ) **2 ) + ( (H ( XU 1 ) - HXMM > * *2 ) ) )*SQRT ( 1 .0 + (DH( 
2 XU 1 )#*?) ) 

2+GW2*HAN2 (G*SQRT( ( ( XU2-XMM) **2 ) + ( (HO.U2 J-HXMM )* *2 )) J*SQRT(1.0+( DH ( 
2 XU 2 )**2) ) 

2+GW3*HAN2<G*SQRT( ( ( XU3-XMM I ** 2 ) + ( (H ( XU3 J- HXMM ) * *2 ) ) J*SQRT( 1.0+(DH( 
2 XU 3 ) **2 ) ) 

2+ GW4 *HAN2 ( G* SQR T ( ( l XU4-XMM ) **2 J + l (H ( XU4 J-HXMM ) **2 ) ) ) *SQRT { 1 * 0 + ( DH ( 
2 XU1>**2)) 

2+GW^*HAN2<G*SQRT( I ( XU5-XMM ) **2 > + ( (HL.U5 J-HXMM >* *2 > > ) *SQR T ( 1 . 0 + ( DH ( 
2 XUS)**2> ) ) 

3361 CONTINUE 

THIS COMPLETES THE FILLlN OF THE MATRIX 
■ NCNSYMMETRIC CRQUT 
FIRST COLUMN OK 
T$0 GET FIRST ROW 
DO 10 J=2 ,N 
10 C(l|J)-C(l«J)/C(l,l) 

C NOW WORK ON ROW AND COLUMN SET < 

DO 11 K=2,N 

KM0=K-1 

KPO=K+l 
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C TO GET DIAGONAL ELEMENT 

S=CMPLX(0. 0,0.0) 

DO 12 l K=1 » KMO 

12 S=S+C(K, IK)*C< IK,K) 

CtK f K ) — C ( Kt K ) -S 

C ' TO GET ELEMENTS IN COLUMN K BELOW ROW K 

IF (KPO.GT.N) GO TO 17 
DO 13 I RO W=KPO» N 
S=CMPLX(0. 0,0.0) 

DO 14 J J= 1 * KMO 
14 S=S+C( IR0W,JJ)*C( JJ ,K) 

13 Cl IROW,K)=C( IR0W,K)-S 

C TO GET ELEMENTS IN ROW K TO THE RIGHT OF COLUMN K 

DC 15 ICOL=KPO,N 
S=CMPLX(0. 0,0.0) 

DO 16 JR= L , KMO 

16 S=S+C(K,JR)*C(JR,ICOL> 

15 CtK, ICOL ) = l C ( K» I COL )-S)/C(K»K) 

17 CONTINUE 
11 CONTINUE 

WRITE (6,1222) N,WE 
1222 FORMAT ( 3H N=,I3,4H WE=,E15.S> 

THI =3. 141 592 7*60. 0/1 80.0 
THXXD=TH 1*1 80.0/3. 1415927 
WRITE (6,9333) THXXO 
9333 FORMAT ( 9H INC ANG=,E15.8) 

THI IS THE ANGLE OF INCIDENCE MEASURED FROM THE HORIZONTAL 
I.E, THE POSITIVE X AXIS 
STH=SIN(TH1 ) 

CTH-COS (THI ) 

THIS FINDS THE INCIDENT FIELD ION THE NJTH SEGMENT 
DO 455 NJ=1,N 
XG=XM{ NJ) 

F(NJ)=CEXP(CMPLX(0.0,G*( ( XG*CTH H- (H ( XG ) *STH » ) )) 


TAPERD ILLUMINATION 


IFIXG.LE. ( (WE*1 .0 )-EP) ) F ( N J ) =CMP LX ( 0 . 0, 0 .0 ) 
l F ( XG .GE . ( EP- ( 1 »0*WE ) ) ) F l N J ) =C MP LX (0 . 0 , C . 0 ) 

IF( (XG.GT.t ( 1.0*WE )-EP ) ) .AND.( XG. LC .( ( 2.0*WE )-EP J ) ) 

2 F(NJ) = FINJ)*(0.5M0.5*S1NI (G/2.0 )*(XG - ( l 1 . 5*WE >-EP ) »)) > 

IPUXG .GE. ( EP-( 2 .0*WE ) ) ) .AND. (XG . LT . I EP-( 1 . 0*WE ) ) ) ) 

2 F{NJ)=F(NJ)*( 0.5-(Q.5*S IN( (G/2.0) * ( XG - ( EP- 1 1 . 5*WE ) ) ) ) ) ) 

ADSF=CABS(F(NJ ) ) 

WR I TE ( 6 , 83 ) NJ.ABSF 

83 FORMAT ( * INC FIELD AT XM ( • , 14, • ) , El 5. 8 I 
455 CONTINUE 

THIS BEGINS THE BACK SUBSTUTION 
CONVERSION CF SOURCE SIDE 
F(l)=F(l)/C(l,l) 

DO 90 I J = 2 » N 
S=CMPLX(C. 0,0.0) 

I JMO= I J-l 
DO 91 I K— l , I JMO 
91 S=S+C( IJ, IK)*F( IK) 

90 F< I J) = ( F(IJ )-S) /C( I J, IJ) 

C NOW FOR FINAL BACK SUBSTITUTION 
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NM0=N-1 

OG 160 L=1 » NMO 

K=N-L 

KP0=K+1 

S=CMPLX<0. 0,0.0) 

DO 175 JO=KPG,N 
175 S=S+C( K,JO)*F( JO) 

160 F(K)=F(K)-S 

DO 425 KCURR=1,N 
ABF=CABS( F( KCURR) ) 

ANGF = 180.0*ATAN2UIHAG(F (KCURR) ) , REAL l FI KCURR ) ) ) /3. 1415927 
425 WRITE(6,553) KCURR, ABF, ANGF 
553 FORMAT ( ' F ( 14 , ' ) = « ,E 15 .8 , • AT. ANGLE* ,E 15 .8 ) 

C THIS ENDS THE BACK SUBSTITUTIONS 

DC 439 KURR=1,N 
IN0=KURR-1 

Y(1)=CABSCF(KURR) ) *4 . 0* WE/ l 6 . 28 3 18*377.0 ) 

XCRD=FLQAT( KURR ) 

439 CALL PL0T(X0RD,Y,1 , IND, 0.0200,0. 0 ) 

DU 440 KURR=1,N 
IND=KURR-1 

Y ( 1 ) =1 80.0*ATAN2 C A I MAG (F{ KUKR) ) , REAL ( FI KURR ) ) )/3. 1415927 
XORD=FLOAT ( KURR ) , 

440 CALL PLOT(XGRD,Y, 1, IND, 180.0,-180.0) 

DC 317 JNX=1 ,360 

TF=0. 87266463 E-02*FL0AT( JNX) 

T=CMPLX(0. 0,0.0) 

DO 310 1=1, N 
Xf«=XM(I) 

310 T=T+ (I F( I )*CEXP(CMPLXI0.C,G*1 <XN*COS(TH) )+( H t XN ) *S IN CTH )))))) ) 
T=T*DC*SQRT (WE) *CMPLX ( -0 .707107 , -0. 707107 ) /3. 1415927 
C)A=C ABS C T ) 

DB=20.0*ALOG10 ( CM) 

CMC ( JNX ) =CM 

CANG=57.296*ATAN2 ( AIMAG(T) ,REAL(T)> 

THD=TH*57.296 

317 WRITE (6,312) Cty, CAHG , THD , DB 

312 FORMAT (18H RELATIVE E F I ELD= , E 15 .8 ,7H ANGLE= »E 15.8, 

2 23H ANGLE FROM HORI ZONT AL= ,E 15. 8 ,6H DB= ,E15.8) 

DO 441 I ES = 1 ,360 

IND=IES-1 

Y C 1 ) =CMC( IES) 

THS=FLOAT( IES) /2.0 
441 CALL PL0T(THS,Y,1, IND, 50. 0,0.0) 

STOP 

END 

FUNCTION H { X ) 

C THIS DEFINES THE SURFACE 

COMMON /PIG/ AONE, CQNE, PONE, ATWO, CTWO, PTWO, N 
H=AONE*SIN(CDNE*X+PONE)+ATWO*SIN( CTWO*X+PTWO) 

RETURN 

END 

FUNCTION HAN21X) 

C I DO THIS TO AVOID RETYPING THE WHOuE GAUSS INT. PART 

COMPLEX HAN2 
COMPLEX AHAN20 
HAN2=AHAN20(X) 

RETURN 

END 


FUNCTION DH ( X ) 

C DHl X) IS THE DERIV. OF H(X) 

COMMON /PIG/ AONE, CONE, PONE, ATVJO, CTWO, PTWO, N 

DH=AGNE*CONE*COSICONE*X+PDNE)+ATWQ*CTWO*COS(CTWO*X+PTWO> 

RETURN 

END 
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C THIS IS A METHOD OF MOMENTS SOLUTION FOR BISTATIC SCATT TM CASE 

C USING T WO POINT INTERPOLATION 

C GAUSSIAN INTEGRATION IS USED TO CALCULATE THE MATRIX ELEMENTS 

C NSUB SEGMENTS HAVE N MIDPOINTS 

C NSUB IS THE SUBSCRIPT WHICH COUNTS THE END POINTS 

C N IS THE SUBSCRIPT WHICH COUNTS THE MIDPOINTS 

C WATCH MAX SLOPE. SD THAT THE X INCREMENTS ARE SMALL ENOUGH 

C THE SURFACE UNpER CONSIDERATION LIES BETWEEN -EP AND +EP 

COMPLEX SNN , SST 
COMPLEX S,CO 
COMPLEX FSS 
COMPLEX FINC120) ,S TS 

COMMON /PIG/ AONE*CQPiE, PONE, ATV/0,CTVf0, PTWO ,N 
COMMON /HOG/ XK( <*00 ) , X { 400 1 , GA , G t DC 

COMMON /G ASSN/ GU1.GU2 ,GU3 ,GU4,GU5 ,GW1 ,6W2 ,GW3 ,GW4,GW5 
COMPLEX C ( 1 50, 1 50 ) 

COMPLEX F<400) , FP( 400 ) , SS , T , CT ES T 
COMPLEX FIN 
COMPLEX HAN2 

DIMENSION A&ES1360) ,YUO) 

C Wg IS THE ELECTRICAL WAVELENGTH 
WE=25.0 

C THE FOLLOWING CONSTANTS DESCRIBE THE SURFACE 

AONE= 15.0 

CONE =2. 0^3. 141 5927/200. 0 
PDNE=0. 0 
ATWO=O,0 
CTWO=0 .0 

ptw<j=o.O 

dc=we/io.o 

OX=DC/iOOO.O 

JDC2=DC/2.0 

DPI F=0 .7853982 

G=6. 2831 853 /WE 

Eb=2. 71828 

GA=G*CC/ (2.0*6E ) 

C EP I S THE END POINT 

EP=200 .0 

C CONSTANTS FOR GAUSSIAN INTEGRATION 5 TH ORDER 

GUI-— 0. 9061798 
GU2=-0. 53846921 
GU3=0 *0 
GU4=-GU2 
GU5=— GUI 
GW1=0. 2369268 
GW 5=0. 2369268 

GU4=0. 47862867 
GW2=0 .47862867 

6W3=0. 5683888 

C CONSTANTS FOR GAUSSIAN INTEGRATION 5 TH 'ORDER 

C THE FOLLOWING BREAKS THE SURFACE INTO SEGMENTS DC CENTIMETERS LONG 

C BY LINE INTEGRATION USING STEPS QF LENGTH DX FOR THE INTEGRATION 

NSUB=1 
X(NSUB)=-EP 

1002 AL=0.000 

t\=X (NSUB) 

1001 R=R+DX 
ALO=AL 

AL=AL+ (DX*SQRT 1 1 .0 + ( DH(R )**2 ) ) ) 

I F ( ( (DC2-AL) .Lg.0.0 ) .AND. ( ( PC2~AL0 ) .GT .0.0) ) XM(NSUB)=R 
IF(AL.LT.DC)GO TO 1001 
NSU8=N5UB+1 
X ( N S UB ) =R 

IF (R.LT.gP) GO TO 1002 
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N=N$UB-i 

• DO 1004 J=1,MSU5 

IP (J.EQ.NSUB) XfA(NSUB)=0.0 
XXX=X( J) 

XMD= XM(J) 

1004 WHITE (6 » 1003 ) XXXtX^D»J 

1003 FORMAT (6H X{ J ) = ,t 1 5 .8 ,9H XM ( J ) = , E 15 3H J=,I3) 

C THIS ENDS THE SURKACb SUBDIVISION 
C THIS INSURES THAT N IS ODD 
KK=0 

5733 KK=KK+1 

IP ( (2*<K-1).EQ.N) GOTO 5731 
IF ( 2*KK.EQ. N) GOTO 5732 
GD TO 5733 
5732 N=N-1 
5731 CONTINUE' 

WRITE (5,3728) Al,KK 

37 28 FORMAT* 1 ‘ , 1 CORRECTED VALUf OF N= * t 1 4 , • KK= • , 14, • 2*KK-1 =N » ) 
NM0=N-1 
NM3=N-3 

DIMENSION OF pIMC,F IS N 
MATR IX F ILL IN 
DO BY COLUMNS 
FOR FI RSI COLUMN 
00 3561 1=1, KK 

3661 CII ,1) =C0 (2*1-1 ,1 K( CO (2*1-1, 2)/ 2.0) 

FOR LAST COLUMN 
DO 3676 1=1 , KK 

3 678 C( 1,KK ) = (CO (2*1- 1,2 *KK-2)/ 2.0) +C0 (2*1-1, 2*KK-1). 

FDR MIDDLE COLUMNS 
DO 56 1 = 1, KK 
11=2*1-1 
KKf11=KK-l 
DO 56 J=2,KKMV 
JJ=2*J-1 

C( I ,J)=(CO(II, JJ-l)/2.0)+C0(I I,JJ ) + { C0( I I, JJ+U/2.0) 

56 CONTINUE 

THIS COMPLETES THE FILL IN OF THE MATRIX 
NONSYMMETRIC CR OUT 
FIRST COLLOM OK 
TO GET THE FIRST ROW 
DO 10 J=2,KK 

10 C( t,J)=C(l,J)/C<ltL) 

C NOW WORK ON ROW AND -COLUMN SET K 

DO 11 K=2,KK , 

KMO=K- 1 
KR0=K+1 

C TO GET DIAGONAL ELEMENT 

S=CMPLX(0. 0,0.0) 

DO \2. I K = 1 , KMO 

12 S=SfC(K,IK)*C( 1K,K) 

C(K,K}=C(K,K}-5 

c to get Elements in column k below row k 

IFIKP0.GT.KK) GO TO 17 
DO 13 IRQW=KP0,KK 
$=CHPLX( 0,0, 0.0) 

DO } 4 JJ=l,KtfQ 

14 S = S + C( IR(JW,JJ)*C(JJ,K) 

13 C( IROW»K!=C( IRQW,K)-S 

C TO GET ELEMENTS IN ROW K TO THE RIGHT OF COLUMN K 

DO 15 I COL=KPO» KK 
S=CMPLX( 0.0,0.01 
DO 16 JR = 1,KM0 

16 $=S+C( K, JR)*C( JR, I COL ) 

15 C ( K , ICOL )= ( C(K,ICOL)-S)/C( K,K) 

17 CONTINUE 

11 CONTINUE 
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WRITE (6,1222) KK, WE 

1222 F(JPM AT ( 1 KK=‘,14,» W£=',E15.8) - 

TH =3. 141 592 7*60- 0/180.0 
THDEG=57.2957S*TH 
WRITE (6,9333) THOEG 
9333 FORMAT { .9H - INC ANG=,E15.8) 

C TH IS THE ANGLE UF INCIDENCE FROM THE HORIZONTAL 

STH=SIMTH) 

CTH=CDS( TH) 

C THIS FINDS THE INCIDENT FIELD ION THE NJTH SEGMENT 
DO 455 NJ=1 , KK 
XG=XM( 2*NJ-1 ) 

FP(NJ)=CEXP(CNPLX( 0.0,G*< (XG*CTH)+(HKXG)*STH) ) ) ) 

IF (XG.LE. ( (WE*1 *0 ) -EP ) ) FP l N J ) =C MPL X ( 0. 0 , 0.0 ) 

I F ( XG. GT . ( EP-1 .0 *W £ ) ) FP(NJ ) =CMP LX (0. 0,0.0) 

IF { ( X6.GT. { < 1.0*WE)-EP) ) .AND.( XG.LE. ( (2.''*WE)-EP) ) ) 

2 FP(NJ)=FP(NJ )*(0,5 + (0.5*SIN( (G/2.0 )*(XG - ( ( 1. 5*WE )-EP ) )) ) ) 
IF( ( XQ.GE. <EP-(2.0*WE) ) ) .AND. (XG.LT • (EP-( 1 .C*WE) ) ) ) 

2 FP(NJ )=FP( NJ ) * l 0 . 5-(0. 5*SIN( ( G / 2. 0 ) * ( XG- ( EP- ( 1 . 5*WE ) ) > ) ) ) 
455 CONTINUE 

WR I TE ( 6,9410) (NJ,FP(NJ) ,NJ=1,KK) 

9410 FORMAT ( 1 •, 'INCIDENT FI ELD F INC( ' , 14, ' ) = » , 2E1 5.8 ) 

THIS BEGINS THE BACK. SUBSTUT ION 
CONVERSION OF SOURCE SIDE 
FP(1)=FP (1 )/C( 1, 1) 

DO 90 I J =2 , K K 
S=CMPLX(0#0»0#0) 

IJMOIJ-I 
. DO 91 IK=1 ,1 JMO 
91 S=S+C( 1 J , IK ) *FP( IK) 

90 FP( I J)=(FP(IJ )-S )/C( IJ, IJ ) 

C NOW FOR FINAL BACK SUBSTITUTION 

NM0=KK-1 
DO 160 L=l,NtfO 
K=KK— L 
KP0=K+1 

S=CMPLX (0.0 ,0.0 ) 

DO 175 JQ=KPO,KK 
175 S=S+C (K, JO )*FP ( JO) 

160 FPtKl =FPIK>-S 
KKM 1= KK- 1 

C TO RECONSTRUCT THE CURRENTS 

DO 47 IRA=1 , KKM1 

47 F( 2* IRA ) = ( FP ( IR A ) +FP ( IRA+1M/2.0 
DO 48 1RA=1,KK 

48 F(2*IRA-1)=FP( IRA) 

WRITE ( 6,4970)( ( J,FP( J) ) ,J = 1,KK) 

497 0 FORM AT ( ’ » , » FP ( • , 1 5, • )= » , 2E 15. 8 ) 

WRITE (6,553) (F(K),K=1,N) 

553 FORMAT ( 6H F( K ) = , 2E 1 5. 8) 

C THIS ENDS THE BACK SUBSTITUTIONS 

DO 439 KURR=1 , N 
IND=KURR-1 

Y(1)=CABS(F( KURR ) )*4.0*WE/ (6.28318*377.0) 

XQRD=FLOAT(KURR ) 

439 CALL PLQT(XORO,Y,1, IUD, 0.0200, 0.0) 

DO 440 KURR =1 ,N 
IND-KURR-1 

Y( 1 )=1 80.0 *ATAN2 (AIMAG(F( KURR) ) , RE AL ( F ( KURR ) ) )/ 3. 1415927 



XORD=F LCAT { KURRJ 

440 CALL PLOT < XORl), Y, 1 , IND, 180. 0,-1 80.0) 

DO 317 JNX=1 , 360 
Th=0.0l7^5323«L LGATtJNX 1/2.0 
T = CM PL X( 0.0, 0.0) 

0C1 310 1=1 ,N 
XN=XM(I) 

310 T = T + ( ( F( I J *CE X P ( CMP L X( 0.0, G*( (XH*COS( TH >>+( H( XN)* SIN ( TH )))))) ) 
********** this CORRECTS the OUTPUT TO TRUE ELE. FIELD 
T=T*DC*SORT( we) *CMPLX( -0.707 107,-0. 707 107 1/3.1415927 
CM® CABS ( T ) 

DB=?0.0*ALOG10 (CM) 

CANG=57.296#ATAN2(AIMAG( T» ,REAL(TM 

THD=TH*57.296 

ABES ( JtJx ) =CM 

317 WRITE (6,312) C M ,C ANG , THO , DB 

312 FORMAT (18H RELATIVE E F I EL 0= , 1 1 5. 8 , 7H ANGLE® , E 1 5. 8, 

2 23K ANGLE F RDM HOR I ZONT AL= i E 1 5 . 8 , 6H ()B=,E15.8) 

DQ 9500 JC® 1,360 
Y ( 1 1 = ABES ( JC) 

£®F LOAY ( JC ) /2.0 
IN D = JC— 1 

9500 CALL PL0T(E,Y,1 ,IND,5G.0,n.C ) 

STOP 

END 


FUNCTION CO (MR , MC ) 

COMPLEX CO 
COMPLEX A HAN 20 

COMMON/ G ASSN/ GUI, GU2 , GU3 , 0U4 , GU 5 , GW 1 , GW2 , GW 3 , GW 4, G W 5 
COMMON /HOG/ XM( 400 ) ,X (400 ) ,GA»G, DC 
lF(MR.tfE*NC) GO TO 100 
C0=DC*AHAN2O(GA) 

GO TO ZOO 
IX CONTINUE 

XMtf=XM(MR) 

HXMM=H ( XMM) 

E PL = X ( MC ) 

EPl) = X ( MC + 1 ) 

DVDFEP® l EPU-EPL) /2.0 
DVSMEP=( EPU + EPL )/2 .0 * 

XU5 = GU5*DVDFEP + PVSME P 
XU 1 = GU 1 *OVDF E P +DVS M E P 
XU2 = GU2*DVOFEP + OVS ME P 
XU3=GU3*DVDFEP+DVSMEP 
XU4=GU4*DVDFEP+0VSMEP 
CO=DVDF£P* ( 

2 + Gkl*AHAN20( G*SORT( ( ( XUl-XMM)**2 »+ ( (H(XU1 )-HXMM) **Z ) ) ) *SQRT(1.C- 
2 + (DH(XUl )**2 ) ) 

2 + GW2*AHAN2 0 < G*SQRT ( ( (XU2-XMM )**2 ) + ( (H( XU 2 )-HXMM )**2 ) )) *SQRT( l. 0 
2 + ( DH ( XU2)**2) ) 

2+GW3*AHANZ0 ( 6*SQRT ( ( (XU3-XMM)**2 )+( (H( XU 3 ) -HXMM ) ** 2 ) ) >*SQRT( l.C 
2 + ( DH ( XU3 )**2 ) ) 

Z+GW^AHANPOl G*SQRT( I ( XU4-XMM ) ** 2 ) + ( ( H ( XU4 ) -HX MM) **2 ) ) )*SQRT(1.0 
2* ( DH (XU4 )**2 ) ) 

2+GW5*AHAN20( G*SQKT ( ( ( XU5-XMM ) *#2 ) + ( ( H ( XU 5 ) -HXMM ) * *2 ) ) J*SQRT( l.C 
2+ ( DH( XU5)*+2») ) 

200 CUNTINUE 
RETURN 
END 


c 


FUNCTION DM ( X ) 

DH ( X ) IS THE DER IV* OF HIX) 

COMMON /FIG/ AfONE»CONf » PONE» ATWQiCTWOf PTWO fN 
DH=AQNE*CON£*COS(CONE*X+PCINE) + AT HG* CTWQ*CGS ( CTWO*X + PTWO ) 


RETURN 

END 


FUNCTION HIX) 

C THIS DEFINES THE SURFACE 

COMMON /PfQ/ AONE, GONE, PONE, AT WO , CTWO ,PT WO ,N 
H=A0NE*S IN(CONE*X+PDNE)+ATWO*SIN(CTWO*X+PTWG) 
RETURN 
END 
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TE CASE .GAUSSIAN INTEGRATION USED TO F ILL IN MA TR I X ,INTE6RfH. 

NSUB SEGMENT.S HAVE N MIDPOINTS 

NSUB IS THE SUBSCRIPT WHICH COUNTS THE END POINTS 
N IS THE SUBSCRIPT WHICH COUNTS THE MIDPOINTS 

WATCH MAX SLOPE SO THAT THE X INCREMENTS ARE SMALL ENOUGH 
THE ARRAY XM(J> CONTAINS THE X COORDINATES OF THE MIDPOINTS OF THE 
SURFACE SEGMENTS ,X(I),X(I+1) ARE THE LOWER AND UPPER X COORDINATES 
OF THE ENDPOINTS OF THE I*TH SEGMENT 

THE SURFACE UNDER CONSIDERATION LIES BETWEEN -EP AND +EP 
cgmplex SNN.SST 
COMPLEX S,CO 
COMPLEX FSS 

COMMON /G ASSN/ GU L , GU2 , GU3 ,GU4 , GU5 , GW 1 , GW2 , GW 3 , GW4 ,GW5 
COMPLEX FINCI20I , S TS 

COMMON /PIG/ AGNE, CONE, PONE, ATWO,CTWO, PTWO, N 
COMPLEX C { 2 35 , 23 5 I 
CCMMON/HOG/ XM ( 400 > , G, X { 400 > 
common /DOG/ DJC 
CCMPLEX CJC 

COMPLEX F1235) ,SS,T,CTEST 
CCMPLEX FIN 
COMPLEX HAN2 

DIMENSION ABES ( 360 ) « Y ( 10 ) 

C THE FQLLOW I NG CONSTANTS DE5CRI BE THE SURFACE 

A0NE=40 • 0 

C0NE=6. 28318/200.0 
P0N£=0.0 
ATWO=0.0 
C T WO = 0 • 0 
P T WO = 0.0 

C WE IS THE ELECTRICAL WAVELENGTH 

WE = 2 5 .0 

G=6. 2831853 /WE 

0C= WE/10 .0 
DX=DC/1OO0.O 
0C2=DC/2 .0 
EP=200 .0 

ST S = -DC + CM PL X ( 0 . 7071 07 , 0 . 707 1 07 ) / ( 2.0 + SORT (WE) ) 

D JC = CMPLX ( C .C,1.0)*G/4.0 

C CONSTANTS FOR GAUSSIAN INTEGRATION 5 TH ' OROER 

GUL=-C« 5061798 
GU2=-C. 53846931 
GU3=0 .0 
GU4=-GU2 
GU5=~GUL 
GW 1=0 • 2369268 
GW5=0. 2369268 
GW4=0. 47862867 

GW 2=0 .47662867 . . 

GW3=0. 5688888 

C CONSTANTS FOR GAUSSIAN INTEGRATION 5 TH ORDER 

C THE FOLLOWING BREAKS THE SURFACE .INTO SEGMENTS DC CENTIMETERS LONG 

C BY LINE INTEGRATION USING STEPS OF LENGTH OX FOR THE INTEGRATION 

NSUQ=1 
X ( N SUB ) =— E P 
1002 AL=0 • COO 
R=XI NSUB ) 

1001 R=R+DX 
ALO=AL 

AL=AL+( DX+SQRTI 1 .0+ ( DH ( R ) **2 )J ) 

I F ( ( ( 0C2-AL ) . LE • 0 . 0 ) . AND . ( (DC2-AL0 ) . GT . 0 .0 ) I XM(NSU8)=R 
IFIAL.LT.DCIGO TO 1001 
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noon 


WRITE (6, 352 ) AL, NSUB 
352 FORMAT { • AL=»,£15.8,' NSUB=',I4) 

NSUB=NSUB+1 

X(NSUB)=R 

IF (R.LT.EP) GO TO 1002 
N=NSUB-i 

V*RITEU,251) NtNSUB 

251 FORMAT l » N=' r I4,» NSUB=*,I4) 

DO 1004 J= 1 , NSUB 
IF (J.EQ.NSUB) XM( NSUB >"0.0 
XXX=X( J) 

XMD= XM { J ) 

1004 WRITE (6,1003) XXX»XMD,J 

1003 FORMAT (6H X ( J ) = , E 15 . 8 , 9H XM ( J ) = ,E 1 5 . 3 , 3H J= , I 3 ) 

C THIS ENDS THE SURFACE SUBDIVISION 

NMC=N— 1 
NM3=N-3 

C DIMENSION OF F!NC,F IS N 

DPIF=0. 7653982 

C MATRIX FILL IN 

00 3661 IR=1 , N 
DO 3661 IC= 1 , N 

3661 Cl IR, IC )=CQ( IR, I C I 

THIS COMPLETES THE FILLIN OF THE MATRIX 
NCNSYKMETRIC CROUT 
FIRST COLUMN OK 
TO GET THE FIRST ROW 
DO 10 J=2 » N 

10 t(l,J)=CIl,J>/CIl,l) 

C NCW WORK DN ROW AND COLUMN SET K 

DO 11 K=2,N 
KMC=K-1 
KP0=K+1 

C TO GET DIAGONAL ELEMENT 

S-CMPLX( 0.0,0.01 
DO 12 IK=1,(<H0 

12 S-S+C ( K , IK)*C( IKiK) 

C(K K J =C ( K » K ) — S 

C 1 TO GET ELEMENTS IN COLUMN X BELOW ROW K 

IF (KPO.GT.N). GO TO 17 
DO 13 IROW=KPQ,N 
S=CMPLX( 0 .0,0.0) 

DO 14 J J = 1 , KMC 

14 S=S+C (IROW , J J)*C( JJ,K) 

13 ci irow,k)=c( mow,K>-s 

C TO GET ELEMENTS IN ROW K TO THE RIGHT OF COLUMN K 

DO 15 )CDL=KPQ,N 
S=CMPLX( 0.0, 0.0) 

DO 16 JR = 1 , KMO 

16 S=S+C(K, JR)*C( JR, ICOL) 

15 C(K, I COL ) = IC(K, ICOL )-S)/ClK,K) 

17 CONTINUE 

11 CONTINUE 

C THIS CNDS THE MATRIX FACTORIZATION 

WRITE (6,1222) N,WE 

1222 F0RMAT(3H N= , 1 3 , 4H WE = ,E15.8) 

TH 1=60. 0*3. 14159/180.0 
WRITE (6,9333 ) .THI 

9333 F0RMAT(9H INC ANG=,E15.8) 

C THI 1 5 THE ANGLE OF INC. MEAS. FROM THE +VE X AXIS 

STH=SIN(THI) 

CTH=COS( THI ) 

C THIS FINDS THE INCIDENT FIELD ION THE NJTH SEGMENT 

DO 455 N J=1 i N 
XG=XM( NJ> 
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THE SIGN ON THE INCIDENT FIELD HAS BEEN ADJUSTED TO AGREE WITH 
THE INTEGRAL EQUATION 

F(NJ)«CEXP(CMPLXtO.O,GM (XG*CTH) +t HI X(,)*STH) ) ) )*CMPLX (-1 .0,0 .0 ) 


TAPPERED I LLUM I NAT ICN 

IF(XG.LE.( (WE*1.0)-EP) ) F < NJ ) =CMPL X( 0 . 0, 0 .0) 

IFCXG.GE. (EP-( I.O*MEm F ( NJ ) =CMPL X( 0 . 0 , 0 . 0) 

IF ( ( XG.GT. ( ( 1 * 0*WE ) - EP J ) .AND . ( XG .L E. U2.Q*WE )-EP) )) 

2 F(NJ)=F(NJ ) *( 0 • 5+ (0 .5*$ INI ( G/2 . 0 ) *( XG - ( { 1 . 5 *WE) -EP ) ) ) ) ) 

IFKXG .GE . (EP-(2.0*WE) ) ) .AND. I XG .LT. { EP- ( I .0*WE ) ) ) > 

2 FIN J I = F( NJ )*( 0 .5- (0.5*SIN( I G/2.0 )* { XG - I EP- ( 1 . 5*WE J ) ) ) ) ) 

CONTINUE 

WRITEI6, 29481 (NJ, F( NJ) ,NJ=1 ,N) 

FORMAT!' INC FIELD F ( ' , I A, ' ) = • ,2E15.8) 

THIS BEGINS THE BACK SUBSTUTION 
CONVERSION CF SOURCE SIDE 

Ft I)=F(1)/C( 1,1) 

DO 90 I J=2 , N 
S=CMPLX( C. 0,0.0) 

I JMO = I J-l 
DO 91 I K = 1 , I JMO 
91 S-S+C( IJ , I K ) *F ( IK) 

90 F( I J ) = C Ft IJ)-S)/C( IJ, IJ) 

C NOW FOR FINAL BACK SUBSTITUTION 

NMC=N-1 

DO 160 L = 1 , NMO 

K=N-L 

KPC=K+ 1 

S=CMPLX(0. 0,0.0) 

DO 175 JO=KPQ,N 
175 S=S+C(K, JQ)*F( JO) 

160 F(K)=F(K)-S 

C THIS ENDS THE BACK SUBSTITUTIONS 

DO 55-4 1KUR=1,N 
A AF=C ABS ( F ( IKUR) ) 

A#F=57.296*ATAN2 ( A I MAG ( F I IKUR) ), REAL l Ft IKUR) ) ) 

554 WRITE.(6, 553) IKUR, AAF,ANF 

553 FORMAT (' ' , * F ( ' , 14, ' ) = ' , E15.8,' AT ANGLE = ' , E 1 5 . 8 ) 

DO 9553 I RRO= I , N 

ind= I pro- i 

Y ( 1 ) -CAB S I F ( IR.RQ ) ) 

XRF0=FL0AT< IRRO) 

9553 CALL PL 0T ( X RRO , Y , 1 , I ND , 5 . 00 , 0 . 0 ) 

DO 9554 I RR 0 = I , N 

IND=TRRO- 1 

Y(1) : »57.2958*ATAN2(AIMAG(F( IRRO) ) , RE AL ( F ( I RRO ) ) ) 

XRRO'FLDAT ( IRRO) 

9554 CALL PLOT ( XRRO ,Y , 1 , I ND, 180.0 ,-180. C) 

DO 317 JNX= 1 , 369 
THS=0.01745329*FL0AT( JNX)/2.0 
T=CMPLX(0.0,0.0l 

DO 310 1=1, N 
XN=XM( I ) 

THN=1.5707963+ATAN(DH(XN) ) 

310 T = T+ HF ( I )*CEXP(CMPLX (0.0 »G*( (XN*COS(THS) ) + (H(XN)*S IN(THS) ) ) ) ) ) 
2 *COS< THN-THS » ) 

C ********** THIS CORRECTS THE OUTPUT TO TRUE MAG. FIELD 

T=T*STS 
CM=CABS( T ) 

DB=20.0*ALOG10ICM> 

. CANG = 57. 296*ATAN2( AIM AG (T) »R E AL ( T ) ) 

THS0=THS*57.296 

ABES(JNX)=CM 

317 WRITE (6,312) CM ,C ANG, THSDtDB 



312 FORMAT <18H RELATIVE H F IELD = , E 1 5 . 8, 7H ANGLE= , E 1 5 .8, 
2 23H ANGLE FROM HQRI ZONTAL=,E15.B,6H DB=,E15.8) 

DO 9500 JC=1,360 
Y ( 1 ) =ABES I JC ) 


U=FLOAT( JCI/2.0 
IND= JC-1 

9500 CALL PLOT<U,Y,l,IND, 50. 0,0.0) 
STOP 
END 


FUNCTION H ( X ) 

C THIS CEFINES THE SURFACE 

COMMON /PIG/ AQNE, CONE ,PONE, ATWO.CTWO, PTWO ,N 
H=AQNE*SIN( (CONE*X ) +PONE ) +AT WO*S IN ( ( CTWO*X)+PTWO ) 
RETURN 
END 


FUNCTION OH { X ) 

C DH(X) IS THE DERIV. OF H(X) 

COMMON /PIG/ AONEtCDNE , PONE, ATWO,CTWO, PTWO, N 
UH^AONE*CONE*COS ( (CONE*X ) +PONE ) +ATWO*CTWO*CO S ( ( C TWO*X ) +PTWO ) 
RETURN 
END 

FUNCTION CO(MR.MC) 

C THIS GIVES THE OLD MATRIX COEFFICIENTS 

COMPLEX CO 
COMPLEX DJC 

COMMON/GASSN/ GU 1 , GU2 , GU3 , GLK ,GU5 , GW1 , GW2 , GW 3 , GWA ,GW5 
COMHCN/HCG/ X M ( 400 ) , G, X ( 400 ) 

COMMON /DOG/ DJC 
COMPLEX AHAN21 
IF(MR.NE.MC) GO TO 100 
C0=CMPLX(0. 500,0.0) 

GO TO 200 
100 CONTINUE 

XMM=XM ( MR ) 

HXMM=H( XMM) 

EP L = X ( MC ) 

EPU=X ( MC+1 ) 

CVDFEP=( EPU-EPL) /2.0 
DVSMEP=( EPU+EPL ) /2 .0 
XU5 = GU5*DVDFEP + OV$ME P 
XU1=GU1*0VDFEP+DVSMEP 
XU2=GU2*DVDFEP+DVSMEP 
XU3=GU3*DVDFEP+DVSMEP 
XU4 = GU'f*DVDFEP + DVSMEP 
HXUl = H( XUI ) 

HXU2=H( XU2 ) 

HXU3=H( XU3 ) 

HXU4=H{ XUA) 

HXU5=H( XU5 ) 

DHXU1=DH(XU1 ) 

DHXU2=DH( XU2 ) 

DHXU3=DH ( XU3 ) 

0HXU4=DH{XU4) 

DHXU5=DH { XU5 ) 

CO=DVDFEP*( 

2+(GWl*AHAN21 ( G*SQRT ( ( ( XU l-XMM ) ** 2 > +( ( HXU 1-HX MM ) * *2 ) ) ) 

2 *( ( -DHXUl# t XMM— XUI ) ) + ( HXMM-HXU 1 ) ) 

2/SQRTU ( XMM- XU L ) **2 > + (( HXMM-HXU D* *2)) ) 

2+( GW2*AHAN21 IG*SQRT( ( ( XU2-XMM ) **2 ) +( (HXU2“HXMM)**2) ) ) 

2 *( ( -DHXU2*( XMM-XU2 ) )+ (HXMM-HXU 2) ) 
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2/SORT ( { ( XNM-XU2 )**2 }+( ( HXMM-HXU2 ) * *2 ) 1 1 
2+(GW3*AHAN21(G*SQRTI ( (XU3~XMM)**2) +( ( HXU3-HXMM J * *2 ) ) ) 
2 *( ( -DHXU3* ( XMM— XU3 ) ) + ( HXMM-HXU3 ) ) 

2/SQKTU (XMM-XU3)**2)+( ( HXMM-HXU3 1**2 1) ) 

2 + ( GW 4* AH AN 21 ( G*SQRT( ( ( XU4-XMM ) ** 2 ) +( (HXU4—HX MM 1**21 I 1 
2 *( ( — DHXU4* ( XMM-XU41 ) + ( HXMM-HXU4 1 1 

2/SQRT ( ( ( XMM-XU4 1 **2 ) + ( ( HXMM— HXU4 1**2 ) 1 ) 

2 + ( GW5*AHAN2l (G*SQRT<{ ( XU5-XMM 1 **2 1 + ( < HXU5-HXMM) **2 1 > 1 
2 *( l-DHXU5*< XMM-XU5 11 + ( HXMM-HXUS 1 I 

2/SQRT (( ( XMM- XU5 1**2 ) + < < HXMM-HXU5 ) * *2 1 1 ) > 

CO=CO*DJC 
200 CONTINUE 
RETURN 
END 


o o o n o o o 


THIS IS TE CASE USING TWO POINT INTERPOLATION 
THIS PROGRAM USES GAUSSIAN INTEGRATION TO GET MATRIX ELEMENTS 
NSU0 SEGMENTS HAVE N MIDPOINTS 

NSUB IS THE SUBSCRIPT WICH COUNTS THE END POINTS 
N IS THE SUBSCRIPT WHICH COUNTS THE MIDPOINTS 

WATCH MAX SLOPE SO THAT THE X INCREMENTS ARE SMALL ENOUGH 
EP IS THE END PAINT 
COMPLEX SNN,SST 
COMPLEX S,CO 
COMPLEX FSS 

COMMON/ GAS SN/ GUI ,GU2, GU3 ,GU4,GU5 ,GWl ,GW2 , GW3 , GW4 , G W 5 
COMPLEX F I NCI 20) , STS 

COMMON /PIG / AON E, CONE, PUNE, AT WO , CT WO ,PTWO » N 
COMPLEX C( 150,1 50) 

COflMON/HOG/ xM< 400) , G, X( 40«) 

COMMON /DOG/ DJC 
COMPLEX DJC 

COMPLEX PI400) , FP< 400) , SS, T,C TEST 
COMPLEX FIN 
COMPLEX HAN 2 

DIMENSION ABES! 360) , Y( 10) 

C WE IS THE ELECTRICAL WAVELENGTH 

WE=25.C 

G=6. 2831853 /WE 

AONE = 5 • 0 

C0N5=6. 28318/200. C 
PCNE = O.C 
ATWO=0.0 
’ CTWO=O.C 
PTWD=O.C 
DC = f/E/l 5. 0 
DX=DC/ 1CCC.G 
DC2 = DC/ 2 »C 
EP=^00.0 

S T S-DC*CM PL X( -0.70711 ,-0. 7071 11 /{ 2. 0* SQRT ( WE) ) 

DJC=CMPLX (0 .0, 1 .D)*G/4 .0 

C CONSTANTS FOR GAUSSIAN INTEGRATION 5 TH ORDER 

GUl=-0. 9061798 
GU2 =-0.53846931 
GU3=0. 0 
GU4=-GU2 
GU5=-GU1 
GW1=0. 2365268 
GW5=G. 2365268 
GW4=0. 47862 867 
GW2 = 0. 47862867 
GW3 =0.5688888 

C CONSTATS FOR GAUSSIAN INTEGRATION 5 TH ORDER 

C THE FOLLOWING BREADS THE SURFACE'INTO SEGMENTS DC CENTIMETERS LONG 

C 8Y LINE INTEGRATION USING STEPS OF LENGTH DX FOR THE INTEGRATION 

NSUfl=l 
X(NSUB)=-EP 
1002 ALsQ.OOO 
R=X (NSUB) 

1001 R=R+DX 
ALO=AL 

AL=AL+(CX*SQkT< 1 ,0+{DH(R J**2) ) ) 

IF ( ( (DC2-AL ) . LE.O.O) .AND. ( ( DC2-AL0) .GT.O.C ) ) XM(NSUB)=R 
IHAL.LT.OC )G0 TO 1D01 
WR1TE(6,352) AL.NSU0 
352 FORMAT! ' AL=',E15.8,' 


NSUB = ' , 14 ) 


o o o o o o o n o n o o 


NSUB=NSUB+1 
X ( NSU8 } =R 

IF (R.LT.EP) GO TO 1002 
N=NSUB-1 

KRITE(6,251 ) N,NSUB 

251 FORMAT! » N= • , I 4, 1 NSUB=»,I4> 

DO 1004 J = 1 » MSUB 
IF (J.EQ.NSUB) XM( NSUB ) =0 *0 

XXX.= X< J I 
X MD= XM(J) 

1004 WRITE {6,1003) XXX,XMD,J 

1003 FORMAT (6H X ( J ) = , E 1 5 . 8 , 9H XM{ J ) = , E 1 5. 8 , 3H J=,I3) 

THIS ENDS THE SURFACE SUBDIVISION 
THIS INSURES THAT N IS ODD 
KK = 0 
5733 KK=KK+ 1 

IF ( ( 2*KK-1) .EC.N) GO TO 5731 
IF (2*KK.EQ.N) GU TO 5732 
GO TO 5733 
5732 N=N-1 
5731 CONTINUE 

WRITE 16,3728) N, KK 

3728 FORMAT! • • , 'CORRECTED VALUE OF N= • , 1 4 , • KK= f , 1 4 , ' 2* KK-1 = N* ) 
NM0=N-1 
NM3=N-3 

DIMENSION OF F I NC , F IS N 
DPIF=0 .7853982 
MATRIX FILL IN 
DO BY COLUMNS 
FOR FIRST COLUMN 
DO 3661 1=1, KK 

3661 C( I, 1 )=C CM2*I-1, 1 )+(C0( 2* I- l, 21/2.0) 

FOR LAST COLUMN 
DO 3678 1=1, KK 

3678 C( 1 ,KK)=(C0(2*I-l,2*KK-2)/2.C)+C0( 2*1- 1 ,2*KK- 1) 

FOR RIDDLE COLUMNS 
DO 56 1=1 , K K 
I 1 = 2*1- 1 
KKM 1 =KK-1 
00 56 J=2,KKM1 
J J= 2*J- 1 

C( I , J )= (C0( I I , J J-l 1/2.0) +C0( I I, JJ )+lCO(II ,JJ + 1 1/2.0) 

56 CONTINUE 

THIS COMPLETES THE FILLIN OF THE MATRIX 
NONSYMMETRIC CROUT 
FIRST C0LLOM OK 
TWO GET FIRST ROW 
DO 10 J=2 , K K 

10 C( 1 ,J)=C( 1, J)/C( 1 ,1) 

C NOW WORK ON ROW AND COLUMN SET K 

DO 11 K=2,KK 
KMO=K— 1 
KP0=K+1 

C TO GET DIAGONAL ELEMENT 

S=CMPLX( 0.0, 0.0) 

• DO 12 I K= 1 , KMC! 

12 S=S+C( K , I K)*C( IK, K) 

C( K,Ki =C( K,K)-S 

C TO GET ELEMENTS IN COLUfAN K BELOW ROW K 

I F ( KPO.GT.KK) GO TO 17 
DO 13 IROW=KPO,KK 
S=CHPLX(0 .0,0.0 ) 

DO 14 J J = 1 ,KMO 
14 S=S+C(IftOW, JJ)*C( JJ,K) 

13 C( I ROW,K)=C( IROW,K)-S 

C TO GET ELEMENTS IN ROW K TO THE RIGHT OF COLUMN K 
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DO 15 ICQL=KPG,KK 
S=CMPLX(0.t,C.O) 

DO 16 JR= 1 • KMO 

16 S=S+C( K, JR )*C( JR, ICOL ) 

15 Cl K, ICOL)=(C(K, ICOL )-S?/C(K,K) 

17 CONTINUE 
11 CONTINUE 

WRITE (6,12221 KK,WE 

1222 FORMAT! •' • , * KK=',I4,« WE=',E15.8| 

TH=3. 141 592 7*60. 0/180.0 
THDEG=57.29578*TH 
• WRITE ( 6,9333) THpEG 
9333 F0RMATI9H INC ANG=,E15.8) 

C TH IS THE ANCLE OF INCIOENCE FROM THE HORIZONTAL 

STH=SIN( TH) 

C T H=C OS ( T H ) 

C THIS FINDS THE INCIDENT FIELD ION THE NJTH SEGMENT 


TAPERED ILLUMINATION ********* ******** ******** 

DO 455 N J = 1 , K K 
XG=XK(2*NJ-1 ) 

FP(NJ) =CEXP(CMPLX(0 .0 ,G*( (XG*CTH) + ( H( XG l*ST H ) ) ) ) *CMPL X( - 1 . C, 0. 0 ) 
C INCIDENT FIELD HAS BEEN ADJUSTED TO AGREE WITH INTEGRAL EQTN. 

IF(X6.LE.(<tfE*l.O J-EP ) ) FP( NJ) =C|WPLX( 0. 0, 0. 01 
1F< X6.GT . ( f P-1 .4#VfE) ) FP( NJ ) = C«PLX (O.C ,C.O ) 

IF( ( XG.GT. ( ( 1.0*WE)-EP ) ) . AND . ( XG. LE . ( ( 2 . r * WE ) -EP) ) ) 

2 FP ( NJ )=FP( NJ )*(0 .5 + ( 0 .5*SIN( (G/2.0 )*< XG -(II, 5* WE )-EP)> ) ) ) 

I F ( ( XG.GE. {EP-(2.0*Vie> ) ) .AND. (XG.LT . (EP-( 1 ,0*WE) ) > J 
2 FP (NJ )=FP( NJ ) *{ 0. 5-( 0.5*SIN( <G/ 2.0 )*( XG-( t P- (1.5* W£ )))))) 

455 CONTINUE 

WR1 TE ( 6 .941 Q) (NJ,FP( NJ) ,NJ = 1 ,KK) 

9410 FORMAT (' 't'lNCrDFNT FIELD F I NC ( * , I 4 , *) = * , 2 E l 5 . 8 ) 

C THIS BEGINS THE BACK SU&STUTION 

C CONVERSION OF SOURCE SIDE 

F P ( 1 )=FP(1)/Cl 1, 1) 

DO 90 I'J = 2 , KK 
S=CNPLX(0.0,0.9) 

I JMCJs: I J- 1 
UO SI I K=1 , I JMO 
91 S = S + C( IJ , IK )*FPl I K I 

90 FP( l J) = (FP( I J 1 -S )/C( IJ, IJ ) 

C NOW FOR FINAL BACK SUBSTITUTION 

NMU = K K - 1 
DO 160 L=1 ,NMO 
K =KK- L 
KP0=K+1 

S=CMpLX (0 .0 ,0.0 ) 

DO 175 JD=KPQ , KK 
175 S=S+C(K,J0)*FP<JQ ) 

160 F P ( K > =F F ( K ) -S 
KKM 1=KK- 1 

C TO RECONSTRUCT THE CURRENTS 

DU 47 IF A=1 , KKM1 

47 F( 2*IRA)=(FP( 1RA)+FP( IRA+1) )/2.C 
DU 48 IRA=1»KK 

48 F l 2*1 RA-I ) =FP ( IRA) 

WRITE ( 6, 49 70 ) ( < J , F P ( J ) ) ,J = 1 ,KK) 

4970 FORMAT!* • , * FP ( ' , 15, * )= • , 2E15.8 ) 

WRITE (6,5531 (F(K),K=1,N) 

553 FORMAT (6H F(K)=,2£l5.8) 

DO IfiRO=l ,N 

I ND=I RRO- 1 
Y ( 1 ) = CAB S ( F( IRRO ) ) 

XRPQ=FLOAT ( IRRO) 



9553 CALL PL OT ( XRRO , Y , 1 , I N D , 5 .00 ,0 .'J ) 

DO 9554 IRR0=1,N 
INP=IRR0-1 

Y(l)=57.2 958*ATAM2tAIMAG(Ft IRRO)) , RE AL ( Ft IRRO ) ) ) 

XRRQ=FLQAT< IRRO) 

955A CALL. PL QT ( XRRQ , Y , 1 , IND , I 80 .0 , -1 83 . 0 ) 

C THIS ENDS THE BACK. SUBSTITUTIONS 

DO 317 JMX=1,360 
IHS=6.01745329*F LT3AT (JNXl/2.0 
T =C MP L X 1 0 . 0 » 0. 0 ) 

DO 310 1 = 1, N 
XN—XM ( I > 

THN=1.5707963+ATAN(DH( XN)) 

310 T=T + ( ( F ( I)*CEXP(CMPLX10.0,G#< ( XN *CDS( THS ) ) + ( H ( XN) * S IN ( THS) ) ) ) ) ) 
2 * CCS ( THN-THS ) ) 

C #£<:** **<-!** THIS CORRECTS THE OUTPUT TO TRUE ELE. FIELD 

T=T*STS 
Cty=CABSCT) 

08= 20 »0#A LOG 10( CM) 

CANG=57.296*AT4)d2 ( A I MAG( T ) . REAL IT)) 

THSD=THS*57.296 

ABES(JNK)=CM 

317 WRITE (6,312) Cfl»CANG, THSD, C8 

312 FORMAT < 18H RELATIVE E F t E LD= , E 1 5 . 8 , 7 H ANGLE= , E 15 .8 , 

2 2 JH ANGLE FRO^ HORIZONTAL®, E15. 8, 6H D.1=,E15.8) 

DU 9500 JC=1,360 
Y ( 1) = ABES( JC) 

U=FLOAT (JCI/2.0 
INQ=JC-1 

9500 CALL PLUTUJ ,Y ,1 , IND, 50. C # 0.0) 

STOP 

END 


FUNCTION H ( X ) 

C THIS DEFINES THE SURFACE 

COMMON /PIG/ ACNE, CONE, PONE, ATWQ, CTWQ,PTWQ,N 
H=AONE*SI N( (CONE*X) +PONE) + ATW(J*SIN< ( C TWO*X ) +P T WO ) 
RETURN 
ENO 


C 


FUNCTION DH ( X ) 

DH( X) I S THE DERIV. OF H(X) 

COMMON /PIG/ AONE, CONE , PONE , AT rtD, CTWO , PTWO, N 
DH=AONC*CQNE*CaS ( (CONE*X ) +PCNE ) +AT WG*CTwO*COS( ( 
RETURN 
END 


CTWO*X ) +PTWO 


) 


FUNCTICN CO(MR,MC) 

C THIS GIVES THE OLD MATRIX COEFFICIENTS 

COMPLEX CO 
COMPLEX DJC 

C CMHON/G ASSN/ GUL.GU2, GU3 ,GLK,GU5 ,GWl , GW 2 , GW 3 , GWA , G W 5 
COMMON/HOG/ XMl AOC ) ,G, Xt AOO) 

common /dog/ oj c 

COMPLEX AHAN21 
I Ft MR .NE.MC ) GO TO 100 
CO=CMPLX (0. 500,0.0) 

GO TO 200 
100 CONTINUE 

XMH=XM(MR) 

HXHH=H ( XMM) 
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EPL=X(MC ) 

EPU=X (KC+1 ) 

DVDFEP=(EPU-EPL) /2.0 
OVSMEP=( EPU+EPL 1/2.0 
XU5=GU5*DVDFEP+UVSMEP 
XU 1-GU l’i'DVDF E P + UVSME P 
XU2 = GU2*DV0FEP-» GVSM EP 
XU3=GU3*DVDFEP+DVSNJEP 
XU4=GU4*0VDFEP+DVSMEP 
AJ0Hl=AT AM DH ( XU1 1 1 
ATDH2 = ATAN< DH ( XU2 I ) 

ATDH3=ATAN( DH( XU 3 ) ) 

ATDH4=ATAN( Dh( XU4 ) ) 

A TOH5=A TAN< DH ( XUS H 
HXUl =H( XU 1 ) 

HXUZ=H< XU2) 

HXU3 = H( XU 3 J 
NXU4=H(XU4) 

HXU5=H( XU5) 

CU=DVDFEP«( 

2 + GW 1>' AHAN2 1 ( G*SQKT ( ( { XU 1 - XMM ) **2 ) + ( ( H( XU 1 I-HXMM ) **2 I ) ) * SQRT ( 1 • 0+( 
2DH ( XU1 >**2 > >*( <-SIN(,ATDHl >*< XMM-XU1 ) I + <C0S< AT CHI )*{ HXMM-HXU1 ) ) ) 
2/SORT ( ( (XMM- XU l )?*2)+( (HXf|K-HXUl)**Z) » 

2 + G W2 ^ AHAN 2 1 ( G^S QfcT ( ( ( XU2-XKM) **2 ) + C ( H < XU2 ) - HXMM ) * *2 ) ) )*SQRT( l.C+( 
2DH ( XU2)**2) ) * ( (-SIN (A T0H2. 1 ^ « XMH-XU2 ) ) + (C0S( AT DH2 J * ( HXMM-HXU2 ) J ) 

2 /SORT ( ( ( XMM-XU2 )**2 ) + ( ( HXMM-HXU2 ) **2) ) 

2+GW3*AHAN21 ( G<= SQRT { ( (XU3-XMM )**2 >+( ( H( XU3 )-HXMM)**2 ) ) ) *SQRT( i.D+( 
2DH(XU3)**2) )*( (-$IN(AT0H3 )*( XMM-XU3) ) + ( COS ( AT IJH3 ) * ( HXMM-HXU3 ) ) ) 

2/ SORT ( ( (Xtfrt-XU3 >**2 ) + ( ( HXHM-HXU3 ) **2 )) 

2+GW4< r AHA^2 1 ( G* SORT { ( { XU4-X MM ) **2 ) + ( { H ( XU4 ) -HX MM ) **2 ) ) ) *SQR T ( I .0 + ( 
2DH( XU4 )**2))*<(-SIN ( ATDH4) *( XHN-XU4) ) + (COS( ATDH4 ) * ( HXMM-HXU4 ) ) ) 

2/ SORT ( 1 ( XMM- XU 4 ) **2 ) + ( ( HXtfM-hXU4 > ** 2 ) ) 

2+GW5*AHAfl21 ( G* SQft T( ( ( XU5-XMM) ** 2 ) + ( ( H ( X U5 ) -HX MM ) **2 ) ) )*SQRT( 1.0+( 
2 DH( XU5 )**2 )>*( ( -S IN ( ATDH5 ) * ( XMM- XU 5) ) + ( Cf)S( A TDH5 ) * < HXMM-HX U5 > ) ) 

2/SQRT( ( (XMM-XU5 )**2)+l (HXMM-HXU5) **2 I ) ) 

CO=DJC*CO 
200 CONTINUE 
RETURN 
END 


FUNCTION AMAN21IX) 

C THIS IS THE HANKEL FUNCTION UF TYPE 2 AND OF ,ORDER 1 

DOUBLE PRECISION XD,DX , A1,A2»A3,A4,A5, A6,HJ1, Bl,B2,B3,B4,B5tAHJl T ' 
2TDX,A1»A2,A3,A4,A5,A6,T1, T2 , T3 , T4 ,T 5 r T6 ,T 7 ,DS Q X, 86 
CCNPLEX AH AN 2 1 
DX=DBLC ( X > 

IF IX.GT.3.0) GO TO 200 
XU=CX*DX/9.0D+Q0 
Al=-0. 3l76lD-O3+0.11C9D-04*XD 
A2=0 .0044331 90+00+ Al^XD 
A3=-0. 039 542890+60+ A2#XD 
A4=<?.210935 730+ ^0+A3«XD 
A5=-0. 562499 850+OQ+A4*XD 
A6=Q. 50+00+ A5*XD 
HJl-AO^DX 

Bi=“0*C40C976D+00+0. 0027 8730+00^X0 

B2=0. 3123 S51D+0C+B1 $X D 

B3=-1.3164827D+00+B2*X0 

B4=2 * 16 82 7090+00 + B3 *X 0 

B5 = 0. 2212091D+00 + B4^XD 

B0=-0 *636619 80 + 00+0 5^X0 

AHJi = ( B6/ OX )+F.Jl*DLOG (DX/2. 0*0 *63661977 
AHAN21=CMPLX(SNGL(HJ1) ,-SNGL(AHJl) } 

GO TO 30 0 
200 T0X=3* O/OX 

A 1=0* CO 1 136530 + 00-0. 00020033* TDX ' 

A2=-0 *002495 1 1 P+OO + Al *TDX 

A3=. 0001 71050 + 00+ A2*TDX 

A4= 0*0 16 5966 70+00+ A3* TO X 

A5=G.15$D-05+A4*TUX 

A6= 0, 79788456O+00 + A5*TOX 

T 1 = 0.000798240 + 0*0-0. 00029 1660+00* TDX 

T2 =G.O0O74348D+OO+Tl*TpX 

T3=-0.t>0637879D+90+T2*TDX 

T 4 =0* 000 056 500 +00 +T 3 *TDX 

T5=0. 1249.96 120+00+T4*TDX 

T6=-2. 356194490+ )f'+T5*TDX 

T7=DX+T6 

0 SO X = A 6/0 SQ RT ( DX ) 

AHAN21=CMPLX(SNGL(DSQX*DCOS(T7) ) ,-SNGL ( DSQX*DSIN( T7 ) ) ) 

300 CONTINUE 
RETURN 
END 
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FUNCTION AH/\N20(X) 

C THIS IS THE HANKEL FUNCTION OF ORDER 0 AND OF TYPE 2 

DOUBLE PRECISION XSQ » BIO , B8 ♦ B6 , B4 t B2 ♦ CIO , C 8, C6 f C4, C2 ,D5 , D4 , D3 
2D2tDlfE5*E4»E2»El ,EO, XD»DX,FO ,E3t HJ ,DSX 
COMPLEX AHAN20 
DX=DBLErl X ) 

IF (X.GT.3.0) GO TO. 100 

XSQ=DX*DX/C.9D+01 

B10-=-0. 39444D-02+XSQ*0.21D-03 

B8=0. 04444 79D+00+XSQ*BlC 

B6=~0.3163866D+00+XSQ*B8 

B4=1.2656208D+00+XSQ*B6 

B2=-2.2499997D+00+XSQ*B4 

HJ=1 .0D+OO+XSQ*B2 

C 10=0. 4279 16D-02-XSQ*G. 248460-03 

C8=-0.42612I4D-01+XSQ*C10 

C6=C.25300117D+00+XSQ*C8 

C4=-0.74350384D+00+XSQ*C6 

C2=0.60559366D+00+XSQ*C4 

hY=SNGL(0.36746691D+00+0.6366198D*00*HJ*DLOG(DX/2.0)+XSQ*C2) 
AHAN20=CMPLXtSNGL(HJ) ,-HY) 

GO TO 200 
100 YD=3.0/DX 

L5=-0.72805D-03+XD*Q. 14476D-03 
04=0 • 1 37237D-02+D5*XD 
D3=-0.9 512D~C44-D4*XD 
l)2=-0.552 74CD-02+D3*X0 
l)l=-0.77D-06+02*XD 
FC=0.79788456U+00+XD*Dl 
E5=-0.29333D-03+XD*0. 13558D-03 
'f 4=-0 . 54 125D-03 + E5*XD 
£3=9. 2625 73D-02+E4*XD 
£2=-C.3954D-Q4+E3*XD ■ 

E1=-0.4166397D-01+E2*XD 

EC= t— 0 .7853981 60+00* XD<=E1) +DX 

DSX=DSQRT { DX ) 

AHAN20=CMPLX(SNGL{ FO^DCOS (EO)/DSX)»-SNGL(FO#DSIN(EO)/DSX) J 
200 CONTINUE 
RETURN 
END 



SUBROUriNEPLOK X, Y, N , I ND , YMAX ,YMI N ) 

01 MENS I0NMU19) ,YLABEL(6 ), Y { 1<? f, tf /)RK l 1 0) 

DATA MARK ( 1 ) , M A R K { 2 ) , MARK ( 3 ) ,MARKC5) , MARK (6) ,NARK l 7 ) , MARK { 8 ) ♦ 

2M ARK ( 9 ) » MARK l 10 ) ,MARKt A) / 1H*, 1H. , 1H1 , lrlO, 1HN,1HH,1H1 ,1HZ , 1 H- , 1HX/ 
DATA I8LANK,N0PT , IPLUS/1H ,1HS,1H+/ 

IF UNO) I* 1* 11 
1 WR I T £ ( 6 , 3 ) 

3 FORMAT (1HI//25X,48H0RDER IN WHICH PLOT SYMBOLS ARE USED *.IX0NH1Z 
*-//30X, 39HTHE SYMBOL U>) INDICATES OFF-SCALE DATA//) 

DC7 J=9 » 1 1 Q 

7 M( J) =MARK( 10 ) 

NCCUNT=1G 

SCALE=lOO .0/ (YMAX-YM IN ) 

LLL = (- YMIN=i c SCALE) + ll .5 

DO 8 J = 1 , 6 

R=J-1 

8 YLABELl J) =P*20»0/SCALE+YMIN 
WRITE(6 i9) (YLABEU I ), 1=1,6) 

9 F0RMAT16X,1PE9.2,5(1PE20.2) / ) 

GOTO 132 

11 NCOUNT=NCOUNT+1 
D099J=1, 119 
99 Ml J )= I BLANK 

IFlLLL.GE.il .AND.LLL.LE.110 ) HI LLL )=MARKl 10) 

IF ( NCOUNT- 10) 133*132, 133 

132 DC89J=U, 111, 20 

89 M ( J ) = I P LU S 

133 D02C J= 1 , N 

L=lYl J)-Y«IN)*SCALE+0.5 • 

IF(L 114,17,17 

14 IFl L+10 >15,16,16 

15 M(1)=N0PT 
GOTO20 

16 LL=L + 1 1 

M( LL)=MARK( J) 

GOTO 20 

17 IFIL-10 8)18,19, 19 

18 LL=L+11 

Ml LL)=MARK( J) 

GOT 020 

19 Ml 119) =NOPT 
2D CONTINUE 

IF ( NCOUNT- 10 ) 21 , 2 5, 2 1 
21 WR I TE ( 6 ,24 ) (MIJ ), J = l,119) 

24 FORMATl IX, 119A1 ) 

GOT 02 7 

25 WR I TE (6 , 26 ) ( X , l M< J J , J = 9, 119)) 

26 FORMATl IX, F7. 3 ,111A1> 

NCCUNT =0 
27 CONTINUE 
RETURN 
END 


APPENDIX B 


SOLUTION OF SYSTEMS OF SIMULTANEOUS LINEAR EQUATIONS 

Several direct methods exist which find the solution -vector, 

[X], when the system of equations 

(99) [C] [X] = [B] 

is given. The two methods used here were the square root (or 
Cholesky) method for syircnetric systems, and the Crout method for 
non-symmetric systems (Ref. [33]). Both methods take advantage of 
the fact that a non-singular matrix [C] is equivalent to [L][U], where 
[L] is a lower triangular matrix and [U] is an upper triangular matrix. 
So 
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... 0“ 


U 11 

U 12 "• U 1N 


c n c i 2 **• c in 


*■21 

*22 
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u 22 ••• U 2N 


C 21 c 22 * * * * 
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min (i , j ) 



( 101 ). 

7T 

11 

£ ik 

U kj 

since 




( 102 ) 

£ ik E 0 

if k 

> i 

( 103 ) 

C 

7T 

C_i. 

Ill 

o 

if k 

> j- 


2 

In order to specify [L] and [U],N +N unknowns must be determined. 

2 

Since there are only N equations, (values of C. .)> N unknowns may 

* J 

be specified. In the square root method the diagonal elements 
are assumed equal, i.e., 

U ii = for 1 = 1, •••. N 

which gives the N extra conditions; in the Crout method one set 
of diagonals is specified, namely 

(104) u kk = 1 for k = 1 , • • • , N. 

Suppose that [C] has been broken up into [L][U], then 

(105) [L][U][X] = [B] 
whence by defining 

(106) [R] = [U][X] 
there results 
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(107) [L][R] = [B] 

which has the solution 

i -1 

(108) r. = (b i - 7 d 1k 1-1, --.N 

and the sum is omitted, if i equals 1. Once the [R] vector is 
known the system 

(109) [U][X] = [R] 

is solved by 

N 

(110) x-j = (r, - [ U ik X|> )/o • • for 1=1 , — ,N 

where the sum is omitted if i equals N. Wilkinson (Ref. [ 3fl) 

has shown that most of the error in a solution of Eq. (99) by 

triangularization methods comes from the decomposition of [C] into 

[L][U] and not in the double back substitution (Eqs. (108) and (110)). 

The details of the decomposition of [C] into [L][U] will now be 

considered. For Crout factorization the diagonal elements of [U] are 

2 2 

set equal to unity leaving N equations and N unknowns in the set 
of Eqs.(lQl), (102) and (103), which can be solved as follows: 

k-1 

Oil) * ik = c ik ' I, ‘im V for i=k, - • • ,N 
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k-1 


01 2) Ukj = i (Ckj ' Ji " km for j=k+1 ■' “ •' 


(113) & ik = 0 if 1 < 


(114) u kj = 0 if j < k. 


These equations are used in the order: first column of [L], first 

row of [U] ; second column of [L], second row of [U]; third column 
of [L], ect. In a computer solution the elements of [U] and [L] may 
be written over the original matrix [C] as they are generated. Once 
this is done the matrix becomes 


FACTORED 


£ 11 * u 12 * U 1N 


AND 

STORED 


* £ 22 

• • 

&M*l • • • • • & 


J N1 


m 


and the fact that the diagonal elements of [U] are unity is used only 
in the previously described back substitution portion of the solution. 
If [C] is symmetric then [C] can be factored into 

015) [C] = [U] T [U] 

fp 

where [U] is the transpose of [U] . Equation (101) becomes 
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(116) 


min of (i ,j) 

ij-] U ki U kj " c i j" 


The u. . ‘s are found from 

1 sJ 


(ii7) u n ~ J c n 


(118) 

(119) 


and 


u ij = c ij /u n for 


i-1 


u ii = (c n - k l, u ki ) 1/2 for i=2 >--> N 


i-1 r. . 

020) u 1j = (c.. - U|ci u kj )/ Uii for|J:^!:;; N > N 


(121) u. . = 0 if i > j. 

The value of this method lies in the reduction of storage space 
required for a given N. With the usual Crout method N storage 
locations are required, but the square root method requires N(N+l)/2 
storage locations since only the upper triangular portion of [C] need 
be stored and [U] can be found using only the upper triangular part 
of [C]. 

A small trick is required if this saving is to be realized in 
practice, since in FORTRAN IV the use of the dimension statement 

o 

"COMPLEX C(N,N) U would set aside N complex storage locations for 


130 



the elements of [C] even if only the upper triangular part of [C] 
were to be filled in and manipulated. To economize on storage a 
way was found to load the elements of the upper triangular part of 
[C] into a linear array N(N+l)/2 positions long. It was convenient 
to preserve the double subscript notation for the matrix manipulations 
and use a simple formula to access the proper location in the singly 
subscripted linear array. A symmetric matrix [C] is shown in Fig. 43 
with the elements of the linear array S inserted into the corres- 
ponding locations of [C]. The order of the matrix is chosen to be 
6 for this example. 


S 1 

s 2 

' S 3 

CO 

S 5 

S 6 


s 7 

to 

CO 

S 9 

s 10 

s n 


s 12 

S 1 3 

s 14 

S 15 



s 16 

S 17 

s 18 




s 19 

s 20 





S 21 


Fig. 43. --Storing a symmetric matrix in a linear array. 

Element c^ is stored in position s -j , a-^ in c 2> etc. The element 
c. . (i <_ 3) can be accessed in the following way. The rows above 

* vl 

the i-th row contain N(i-l) - ( (i-1 )(i -2)/2) elements and in the i-th 
row there are j - i+1 elements up to and including the one to be 
accessed, hence 
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(122) c., = s(N(i-l) - il£UiL^)_ + j - i + i) 

* J ^ 

=(sN-i-[(l^i-) + N - j]) . 

In the programs the subscript manipulations are performed directly 
in the subscript or accessed by calling a function named ISUB(i,j) 
[Integer Subscript corresponding to i,j]. If, for example, c-j 5 were 
needed in a computation the element s (ISUB(1 ,5) ) is used. Once the 
factorization is completed, the back substitutions are performed. 

Notice that in either the Crout method or the square root 
method there are two distinct steps. The first is factoring the 
matrix and the second is the back substitution. The first step is 
independent of the driving column [B] and hence need be done only once 
for any given matrix [C] so any number of driving columns may be 
considered without re-factoring [C] . 
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